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STRESZCZENIE 

Niniejsza praca prezentuje nowo opracowany model ewolucji hydrodynamicznej, 
ktory w pol^czeniu z modelem statystycznej hadronizacji THERMINATOR sluzy nam do 
opisu zachowania silnie oddzialuj^cej materii wyprodukowanej w relatywistycznych 
zderzeniach ci^zkich jonow. Nasze oryginalne podejscie wykorzystano do wykonania 
dopasowan dla danych pochodz^cych z eksperymentow realizowanych na akcelera- 
torze RHIC (Relativistic Heavy Ion Collider w Brookhaven National Laboratory) 
przy najwyzszej jego energii ^snn = 200 GeV, oraz do sformulowania przewidywan 
teoretycznych dla przyszlych eksperymentow cigzkojonowych przy wyzszych cner- 
giach (a/snn = 5.5 TeV, dla akceleratora LHC, skrot od Large Hadron Collider w 
CERN-ie). 

Nasze wyniki odnosz^ siq do obserwabli jedno- i dwu-cz^tkowych w zakresie 
mi^kkiej fizyki (px < 2 GeV). Opisujemy widma cz^tek w p^dzie poprzecznym, 
wspolczynnik przeplywu eliptycznego V2, oraz promienie HBT dla identycznych pio- 
now (HBT jest skrotem od nazwisk Hanbury-Brown i Twiss). W ramach prac 
nad rozpraw^ skonstruowalismy nowe rownanie stanu dla materii silnie oddzialu- 
jqcej, ktore l^czy model gazu hadronowego z wynikami symulacji QCD na siatkach. 
Calosc programow tworzy platform^ obliczeniow^ w sklad ktorej wchodzi kod hydro- 
dynamiczny pol^czony z modelem statystycznej hadronizacji THERMINATOR (skrot od 
THERMal heavy loN generATOR). 

Stosuj^c standardowy model optyczny Glaubera jako warunck pocz^tkowy dla 
ewolucji hydrodynamicznej, osiagn^lismy bardzo dobry opis danych eksperymen- 
talnych uzyskanych na akceleratorze RHIC. W szczegolnosci osiagn^lismy znacznie 
lepszy, od wczesniej uzyskanych w modelach hydrodynamicznych, opis promieni ko- 
relacyjnych. Dla przyszlych eksperymentow ci^zkojonowych na LHC uzyskalismy 
przewidywania teoretyczne dotycz^ce mi^kkich obserwabli. 

Zaproponowalismy rowniez sposob rozwi^zania tzw. zagadki HBT na RHIC-u. 
Sugerujemy zmodyfikowanie warunkow pocz^tkowych i wprowadznie Gaussowskiego 
profilu g^stosci materii jako warunku pocz^tkowego dla hydro dynamiki. Taka mody- 
fikacja prowadzi do szybszcj formacji poprzecznego przeplywu kolektywnego, co 
warunkuje uzyskanie wyj^tkowo dobrej zgodnosci naszego modelu z danymi ekspery- 
mentalnymi. 

Jako ostatni punkt, wprowadzilismy do naszego modelu proces swobodnego stru- 
mieniowania cz^stek, ktory w pol^czeniu z mcchanizmcm naglej, chociaz opoznioncj 
w czasie termalizacji, tworzy nowe warunki pocz^tkowe dla kodu hydrodynami- 
cznego. Wprowadzenie przedrownowagowej ewolucji pozwolilo nam na opoznienie 
startu fazy hydrodynamicznej. Jest to poz^dany efekt, ktory ma na celu unikni^cie 
zalozenia o bardzo wczesnej termalizacji ukladu, ktore wydaje si§ bardzo trudne 
do uzasadnienia na gruncie mikroskopowym. Podkreslmy, iz wl^czenie swobodnego 
strumieniowania cz^stek nie zmienia wysokiej zgodnosci uzyskanych wynikow mode- 
lowych z danymi doswiadczalnymi. 
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/ know that this defies the law of gravity, 
but, you see, I never studied law. 

- Bugs Bunny 



Physics isn't a religion. 
If it were, we 'd have a much easier time raising money. 

- Leon Lederman 



— to my Parents — 
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Chapter 1 
Introduction 



Nowadays, the relativistic hydrodynamics is regarded as the best theoretical frame- 
work for description of the spacetime evolution of strongly interacting matter pro- 
duced in ultra-relativistic heavy-ion coUisions PElEHHElinilllElinillOllIIlliailHl 
[lllIlSlliniinilTHlIiniEniEniiaiEaEalE^ recent review see [25] . In particu- 

lar, the soft hadronic one-particle data describing the transverse-momentum spectra 
and the elliptic flow coefficient f2, collected in the RHIC experiments (Relativistic 
Heavy- Ion Collider at the Brookhaven National Laboratory), have been successfully 
explained in various approaches based on the perfect-fluid hydrodynamics. In fact, 
the explanation of the large value of the elliptic flow by perfect hydrodynamics is 
regarded as the evidence of early thermalization and suggests that the quark-gluon 



On the other hand, the approaches based on the hydrodynamics cannot repro- 
duce the two-particle observables such as the pion correlation radii. The latter are 
commonly called the HBT radii - after Hanbury-Brown and Twiss who in 1950s 
showed that it was possible to determine the angular sizes of astronomical radio 
sources and stars from the correlations of signal intensities, rather than amplitudes. 
The difficulty of the consistent description of the one- and two-particle observables 
picked up the name "HBT puzzle". The HBT puzzle and the problem of the mi- 
croscopic explanation of very fast thermalization of the matter produced at RHIC 
represent two issues that challenge the hydrodynamic picture. 

In this Thesis we present our recently developed hydrodynamic model and use it 
to describe the RHIC data. We address both the one- and two-particle observables: 
the transverse momentum spectra, the elliptic flow coefficient f2, and the pion HBT 
radii. We suggest how the HBT puzzle as well as the early thermalization problem 
may be solved. Similarly to other approaches, our model is based on the perfect 
fluid hydrodynamics and includes the symmetry against Lorentz boosts along the 
beam axis, the so called boost-invariance. This restriction means that our results 
may be applied only to the central regions of relativistic heavy-ion collisions. On 
the other hand, our framework differs from other approaches in several important 
aspects, in particular, in the use of a different equation of state, modiflcation of the 
initial conditions, two-body method of the calculation of the correlation functions. 



plasma created at RHIC is a strongly interacting system 
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and different treatment of tlie final fiadronic stage. Tlie main acliievements of tlie 
Tfiesis are tlie following: 

1. The construction of the realistic equation of state for strongly interacting mat- 
ter which interpolates between the hadron gas model and the results of the QCD 
lattice simulations. This equation of state describes the crossover phase transi- 
tion, i.e., the transition where thermodynamic variables such as energy density 
or entropy density change very rapidly in the narrow range of the tempera- 
ture, however, no real discontinuities in the behavior of the thermodynamic 
variables are present. Noticeably, our equation of state has no pronounced 
soft point where the sound velocity is very small and possibly drops to zero. 
It is known that the presence of such a soft point leads to the ratio of the 
HBT radii i?out/-Rside that is larger than unity - an effect which has not been 
confirmed by the experimental data. 

2. The boost-invariant hydrodynamic equations for baryon free matter have been 
rewritten in the very concise form which reduces the number of the independent 
equations to two. This may be done in the formal way if the range of the vari- 
able r (the distance from the collision axis) is extended to the negative values. 
The applied procedure is a direct generalization of the formalism introduced 
earlier by Baym et al. in the studies of cylindrically symmetric systems with 
constant sound velocity. The new form of the hydrodynamic equations allows 
for the simple and natural inclusion of the boundary conditions at the origin 
of the system. 

3. The computational platform has been constructed which combines the hydrody- 
namic code with the statistical hadronization moc?e/ THERMINAT0R[27j. This is 
arranged in such a way that the information about the freeze-out hypersurface 
obtained from the hydrodynamic code is exported and treated as an input 
for THERMINATOR. We emphasize that our equation of state (in the region be- 
low the critical temperature T^ ~ 170 MeV) describes the hadron gas with the 
same set of the hadronic species as that included in THERMINATOR. Hence, there 
is a smooth change between the hydrodynamic and statistical description of 
the produced matter. THERMINATOR is a Monte Carlo program simulating the 
decays of resonances. The Monte-Carlo method allows for the direct compar- 
ison of our model results with the experimental data. In particular, one can 
easily include various experimental cuts. 

4. The successful description of the soft hadronic RHIC data has been achieved 
with the standard initial conditions obtained from the optical limit of the Glauber 
model. By this we mean here that the use of the new equation of state helped 
to reduce the discrepancy between the theoretical and experimental ratio of 
the HBT radii -Rout and -Rside- In this version of our calculations we find 
-RoutZ-Rside ~ 1.25, significantly closer to the experimental values than in ear- 
lier hydrodynamic studies. 
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5. Predictions for the future heavy-ion collisions at LHC have been formulated. In 
the studied by us central region we expect that the transition from the RHIC 
energy, ^Jsnn = 200 GeV, to the LHC energy, ^/sNN = 5.5 TeV, results 
essentially in a higher initial central temperature Ti used as the input for the 
hydrodynamic calculation. Thus, one can make predictions for the collisions at 
the LHC energies using a set of values for T[ which are higher than those used 
at RHIC. At RHIC we found Tj = 320 MeV, hence for LHC we studied the 
cases Ti = 400, 450, and 500 MeV. Our results for LHC indicate a moderate 
increase of the HBT radii and saturation of the pion elliptic flow (as compared 
to the RHIC experiments). 

6. The solution of the RHIC HBT puzzle has been proposed which suggests the 
use of the modified Gaussian-type initial conditions. We find that the choice of 
the initial condition in the form of a two-dimensional Gaussian profile for the 
transverse energy leads to a complete and consistent description of soft observ- 
ables measured at RHIC. The transverse-momentum spectra, the elhptic-flow, 
and the HBT correlation radii, including the ratio i?out/-Rside are very well 
described. 

7. The processes of the free streaming of partons followed by the sudden equili- 
bration were incorporated in the model. Those two processes deliver modified 
initial conditions for the hydrodynamics. In particular, the inclusion of the 
free-streaming stage allows for the delayed start of the hydrodynamic evo- 
lution, which is a desirable effect in the context of the early thermalization 
problem. 

The Thesis is organized as follows. In Chapter 2 we describe the construction 
of our equation of state. In Chapter 3 we present the hydrodynamic equations and 
transform them to the form used in the numerical calculations. The initial conditions 
and the freeze-out prescription are introduced in Chapters 4 and 5, respectively. The 
fits to the RHIC data obtained with the standard initial conditions are presented 
in Chapter 6. The predictions for the LHC are given in Chapter 7. The solution of 
the RHIC HBT puzzle with modified Gaussian initial conditions is presented and 
discussed in Chapter 8. In that Section we also discuss the inclusion of the parton 
free-streaming as the pre-hydrodynamics stage. The Summary and four Appendices 
close the Thesis. 

Wc use everywhere the natural units with c — h — — 1- The signature of the 
metric tensor is (-1 ). 
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Chapter 2 



Thermodynamics of relativistic 
baryon-free matter 

In our approach we concentrate on the description of the mid-rapidity region of 
ultra-relativistic heavy-ion colhsions. Statistical analysis applied to the highest- 
energy RHIC data indicates that the baryon chemical potential hb at the chemical 
freeze-out is of about 25 MeV in this region [281 EH], [301 EH [32] • The predictions of 
the statistical models for LHC give even smaller values, /xb ~ 0.8 MeV [33j- On the 
other hand, the expected temperature is of about 150 - 170 MeV, hence the ratio 
Hb/T is small and in the hydrodynamic equations we can approximately assume that 
the baryon chemical potential vanishes. In this situation, as discussed in Ref. [34j . 
the whole information about the equation of state is encoded in the temperature- 
dependent sound velocity Cs(T). We assume that at low temperatures the sound 
velocity is given by the hadron-gas model with a complete set of hadronic resonances. 
In this case the function Cs(T) approaches zero as T/ttIt:, which is the characteristics 
of the pion gas [771^^ is the pion mass). On the other hand, at high temperatures 
our equation of state coincides with the recent lattice simulations of QCD [35j- 
The thermodynamic properties of the hadron gas and the quark-gluon plasma are 
discussed below in more detail in Sects. 12.11 and 12.21 respectively. 

In the transition region between the hadron gas and the plasma, whose position 
is characterized by the critical temperature Tc, different interpolations between the 
hadron gas result and the lattice result may be considered. In Ref. we showed, 
however, that the most promising equation of state is based on the simplest inter- 
polation between the hadron-gas model and the lattice data, see Fig. 12.11 This is so 
because the sound velocity function which does not exhibit a distinct minimum at 
the critical temperature leads to the relatively short evolution time and this effect 
helps to describe correctly the HBT data. The effects of different forms of the sound 
velocity are discussed in Sect. 12.31 Since the simplest interpolation is the best, most 
of the results presented in this Thesis are obtained with the sound velocity function 
shown in Fig 12.11 

The knowledge of the function Cs(T) allows us to determine all other thermo- 
dynamic properties of our system. This is achieved with the help of the following 
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Figure 2.1: Temperature dependence of the square of the sound velocity at zero baryon 
density. The plot shows the result of the lattice simulations of QCD [35] (solid line) 
and the result obtained with the ideal hadron-gas model discussed in Sect. 12.11 (dashed 
line). A piece of the thick solid line describes the simplest interpolation between the two 
calculations. The critical temperature equals 170 MeV. It is defined as the place where 
the sudden change of the thermodynamic variables occurs, see Fig. 12.21 

thermodynamic identities 

dp 

e + P = Ts, de = Tds, dP = sdT, cl = — , (2.1) 

where e is the energy density, P is the pressure, T is the temperature, and s is the 
entropy density. In Fig. I2.2l we display the entropy and energy densities as functions 
of T, and the pressure and sound velocity as functions of the energy density. These 
quantities follow directly from the assumed form of the function Cs(T). shown in Fig. 

We note that other equations of state and their impact on the physical ob- 
servables were recently studied in Ref. |27]. The result of that work was that the 
transverse-momentum spectra are quite insensitive to the assumed form of the equa- 
tion of state. On the other hand, a noticeable dependence of the elliptic flow on 
the equation of state was observed. This dependence favored the strong first order 
phase transition. In Ref. the effects concerning the HBT were not studied. 
Our recent work indicates that the constraints from the particle interferometry data 
exclude the strong first order phase transition since it leads to the unrealistically 
long evolution times. 

We stress that by using the lattice results we take into account the non - per- 
turbative aspects of the plasma behavior, which may be regarded as the effective 
inclusion of the strongly-interacting quark-gluon plasma; large deviations from the 
ideal-gas behavior directly indicate the non-negligible interactions present in the 
plasma. In particular, Cg is significantly below the ideal-gas value of 1/3 also at 
temperatures way above Tc. Note that in agreement with the present knowledge, no 
real phase transition is present in the system, but a smooth cross-over, therefore Cg 
does not drop to zero at T^ but remains a smooth function. 
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Figure 2.2: The two left panels: the entropy and energy densities, scaled by and T^, 
respectively, shown as functions of the temperature. The two right panels: the pressure and 
sound velocity shown as functions of the energy density. The presented thermodynamic 
functions follow directly from the temperature-dependent sound velocity shown in Fig. 12.11 
We observe a sudden but smooth change of s/T^ and e/T^ at T ~ Tc. The vertical line 
indicates in all cases the critical energy density corresponding to Tc = 170 MeV. With our 
equation of state one finds Ec = 0.6GeV/fm^. 



2.1 Hadron gas 



The hadron-gas model is based on the assumption that all known hadrons (including 
hadronic resonances) form a multicomponent perfect gas. In this case, well known 
formulas for the thermodynamic variables of the relativistic perfect gases 
may be applied. Of course, the hadron-gas model is applicable in the temperature 
region below the critical temperature Tc = 170 MeV. With increasing temperature 
the density of hadrons becomes so large that they start overlapping and the very 
idea of hadrons breaks down. We note that only in the case of pions the quantum 
statistics is necessary. For other particles, which are much heavier and less abundant, 
the Boltzmann classical limit is sufficient 



20 



CHAPTER 2. RELATIVISTIC THERMODYNAMICS 



2.1.1 Pion gas 

Pions are the most abundant particles produced in heavy-ion colhsions. Because of 
their large multiphcity the statistical quantum effects cannot be neglected and one 
should use the Bose-Einstein distributions to describe the pion spectra. Additionally, 
since pions do not carry the baryon number and strangeness, their chemical potential 
may be assumed to be zero (possible values of the isospin chemical potential are 
very small, fii^ < 1 MeV, and are usually neglected). Therefore, in describing the 
thermodynamic properties of the pion gas (ttG) we use the relations for a massive 
boson gas with zero chemical potential, which are worked out in Appendix |X1 namely 

1 OO 

P'°(T) = (2.3) 

ft=l 
1 oo 

efiT) = ^,Tm^J2^2[^TK,{^n)+m.nK,{^n)], (2.4) 



K=l 



where the index i specifies the isospin (distinguishes between vr^, 7r°, 7r+) and rrii is 
the mass of the appropriate pion. The infinite sum over n is the technical way to 
include the Bose-Einstein statistics. In the numerical calculations it is sufficient to 
include only the first four terms - the inverse powers of k, reduce the higher-order 
terms fast enough. 

Knowing the expression for the entropy density of the pion gas, we calculate the 
sound velocity from the last equation in (12. ip . In the approximation that the pion 
masses are equal (mj = m^) we obtain 



oo 

{clY'' (T) = ^^^^ • (2.5) 

K=l 



It is interesting to observe the behavior of the sound velocity in the limit when 
temperature approaches zero. The arguments of the Bessel functions Ki,{z) tend to 
infinity and we can use the power series expansion to find that 



K)^^(T-O) 



^— \oo 

T Z^«=i 


3 

K 26 


T -\- ■ ■ ■ 




1 

K, 2e 


T -\- ■ ■ ■ 



T 



(2.6) 



see Eq. ( ID.SP from Appendix [Dl On the other hand, when temperature is large 
(T ^ m^) we expand the modified Bessel functions according to (]D.7p and get the 
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following expression 



oo 



1 Z]k=i ^ 8 ("r +8(T'^) +■ 

3 V^oo 



24 V T ^ 96 \ T '^J 



1 

3' 



:2.7) 



The last limit illustrates the expected result - at very high temperatures the massive 
pion gas behaves effectively like a massless gas with P = and its sound velocity 
squared equals |. (See Tables and 1X31 where the complete formulas for pressure 
and energy density in various limits are given). 

2.1.2 Classical gas 

Strictly speaking, all particles produced in heavy-ion collisions obey the quantum 
statistics, however for all particles other than pions the effects of quantum statistics 
are numerically negligible. Thus, to calculate the thermodynamic functions of the 
hadrons other than pions we may use the following formulas 



9i 3 (^^ 
7m Aq — 

27r2 ' H T 

-^T^mlKj— 
27r2 * H T 



27r2 



T m- 



rrii 



T 



?>TK2{-^]+m,KA^ 



rrii 



T 



(2.8) 
(2.9) 
(2.10) 



where (^j = 2 Sj + 1 is the degeneration factor which holds the information about the 
spin degeneration of i-th particle. The sound velocity in the classical massive case 
is defined as follows 



K)r(^) 



2T 



m, 



(2.11) 



and has the same asymptotic features as the pion gas, namely it tends to the value | 
for very high temperatures and becomes proportional to T if the temperature tends 
to zero. 



2.1.3 Massive hadron gas 

The massive hadron gas (HG) model is the sum of both massive pion gas and massive 
classical gas of all other hadrons. In our study, the information on the mass and 
spin of individual particles comes from the input file to the SHARE program (41j . 
The table particles, data holds parameters for 371 particles consisting from u, d and 
s quarks. Thus, for all hadrons we may write 



371 



^HG^j.^ = $^(3)f (T) + 5^(.)f^(T), 
(cs^)''^(T) 



dT 



i=A 



T d{sHGy 



(2.12) 
(2.13) 
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and all other thermodynamic quantities (s, P and e) are expressed in the analogous 
way. 



2.2 Quark-Gluon Plasma 

In our approach we use the results of the lattice simulations of QCD at the finite 
temperature presented in Ref. [35], see Fig. 12. 3[ They were obtained for physical 
masses of the light quarks and the strange quark. The pressure data obtained from 
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Figure 2.3: QCD pressure as a function of temperature normalized by T'^ [3S]. Results 
are obtained for two lattice spacings A't = 4 (red) and A't = 6 (blue). 

the QCD lattice calculation in Ref. [35] have been recently parameterized for the 
case of A^t = 6 in Ref. [42], see Fig. 12.41 The parameterization has the form 

1 — — 

P = cT^a{TjT), aig) = ^^e-'<^, (2.14) 

1 + eV 

where the dimensionless fit parameters equal: a = 0.91, b = 0.11, c = 5.21 and 
A = 1.08. Taking the lattice result at the face value, one expects that the sound 
velocity significantly drops down in the region T ^ T^. Similar behavior, with Cs(Tc) 
reaching zero, is expected in the case of the first order phase transition where the 
changes of the energy happen at constant pressure. However, the lattice simulations 
suggest that for three massive quarks with realistic masses we deal with the cross- 
over rather than with the first order phase transition, hence the sound velocity 
remains finite, as is consistently shown in Fig. 12.11 





Figure 2.4: Pressure as a function of temperature normalized to the pressure corre- 
sponding to the Stefan-Boltzmann limit. Points represent results from the QCD lattice 
simulations [35], whereas the line is the data fit from Ref. |42j . 

2.3 Modeling the crossover phase transition 

The exact values of the sound velocity in the region T ^ are poorly known. The 
lattice calculations are not very much reliable for T < and, at the same time, the 
use of the hadron gas model with vacuum parameters becomes unrealistic for large 
densities (temperatures). The authors of Ref. [35] state that in the hadronic phase 
the lattice spacing is larger than 0.3 fm and the lattice artifacts cannot be controlled 
in this region. In this situation, it is practical to consider different interpolations 
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Figure 2.5: Three different forms of the sound- velocity function analyzed in Ref. [36j. 
The solid line describes the interpolation between the lattice and the hadron-gas results 
|34] with a shallow minimum where Cg = 0.14 (case I), the dashed line describes the 
interpolation with a dip where Cg = 0.08 (case II), finally the long-dashed line describes 
the interpolation with a deep minimum where Cg = 0.03 (case III). Note that the case I is 
the approximation of the result shown previously in Fig. 12.11 - see discussion in the text. 
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Figure 2.6: The temperature dependence of the entropy density and energy density, panels 
(a) and (c), as well as the energy density dependence of the pressure and sound velocity, 
panels (b) and (d). One can observe that the deeper is the minimum of the sound velocity 
function, the steeper is the increase of the entropy density and the energy density. 



between the lattice and hadron-gas results and to analyze the physical effects of a 
particular choice of the interpolating function. This type of the study was performed 
in Ref. [36]. Here we shortly discuss the main conclusions of this analysis. 

In Ref. [36] we considered three different sound- velocity functions Cs(T). Below, 
we refer to these three options as to the cases I, II and III, see Fig. 12.51 In the case 
I, we use the sound-velocity function which agrees with the ideal hadron gas model 
of Ref. [M] in the temperature range < T < 0.85 T^, and with the lattice result in 
the temperature range T > 1.15 TJ1|. In the region close to the critical temperature, 
0.85 Tc < T < 1.15 Tc, a simple interpolation between the two results is used. We 
have checked that such a simple interpolation yields directly the entropy density 



^Ref. [34j uses the approximation where the sum over hadronic states is replaced by the integral 
with the mass-density functions worked out in [13]. This leads to small differences between Cs{T) 
used in [31 [35] and [S] 
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consistent with the lattice result. Namely, the use of the thermodynamic relation 



relating the entropy density with the sound velocity for zero baryon chemical poten- 
tial, gives the function s{T) which agrees with the lattice result at high temperatures, 
s{T)/T'^ ^ 12 at T = 1.5 [33]. 

In the cases II and III, the sound-velocity interpolating functions have a distinct 
minimum at T = Tc. Comparing to the case I [with Cs(Tc) = 0.37 and Cg (Tc) = 0.14], 
the value of the sound velocity at T = Tc is reduced by 25 % in the case II [where 
Cs(Tc) = 0.28 and cl{T^) = 0.08], and by 50% in the case III [where Cs(Tc) = 0.19 
and Cg(Tc) = 0.03]. From Eq. (12.151) one concludes that the decrease of the sound 
velocity at T^. leads to the increase of entropy density for high temperatures. Hence, 
in order to have the same value of the entropy density at high temperatures, a 
decrease of the sound velocity function in the region T pa Tc should be compensated 
by its increase in a different temperature range. For our interpolating functions 
in the cases II and III we assume that the values of Cs(T) in the range 0.15 Tc < 
T < 0.85 Tc are slightly higher than in the case I, see Fig. 12. 5[ Such modifications 
may be regarded as the parameterization of the repulsive van der Waals forces in 
the hadron gas. The values of the maxima are chosen in such a way that the 
entropy densities for three considered cases are consistent with the lattice result, 
see the upper left panel of Fig. 12.61 where the functions s(T)/T^ are shown. We 
stress that in the three considered cases the values of Cs(T) in the temperature 
range Tc < T < 1.25 Tc remain significantly below the massless limit l/y/S. Such a 
limiting value is implicitly used in many hydrodynamic codes assuming the equation 
of state of non-interacting massless quarks and gluons for T > Tc, see for example 
the extended 3+1 hydrodynamic model of Ref. |23j . 

The studies of the hydrodynamic spacetime evolution of matter described by 
the equations of state I, II and III were performed in Ref. [SB]. As expected, we 
found that the drop of Cg at T = Tc leads to the prolonged time of the evolution, 
hence leads to the increase of the -Rout/-Rside ratio of the HBT radii. This behavior 
is in contradiction with the observed data which indicate -Rout/-Rside ~ 1- Thus, we 
have decided to exclude the cases II and III from further analysis and to restrict our 
consideration to the case I. 




(2.15) 
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Chapter 3 



Relativistic hydrodynamics of 
perfect fluid 



In this Chapter we present the main ingredients of our hydrodynamic modeL We 
start with the general formulation of the relativistic hydrodynamics in the case of 
vanishing baryon number. Next, we implement the idea of boost-invariance that al- 
lows us to restrict our considerations to the plane z = 0. In the subsequent Sections 
of this Chapter we show that the hydrodynamic equations may be rewritten in the 
form where only two equations are independent (at the expense of the formal exten- 
sion of the variable r = a/x^ + to negative values) and the boundary conditions 
at the origin are automatically fulfilled. Such a form, being the direct generalization 
of the approach introduced in Ref. ^j, turned out to be very convenient in the 
numerical analyses. 

3.1 Hydrodynamic equations 
for baryon free matter 

As the system reaches local thermodynamic equilibrium its further evolution is gov- 
erned by the conservation laws for energy and momentum, which can be expressed 
by the formula 



Here e is the energy density, P is the pressure, g'^'^ is the metric tensor (we use the 
convention where g^^ = +1) and u'^ is the four- velocity. 



In Eq. (13.31) v is the local three-velocity of the fluid and 7 is the Lorentz factor 




(3.1) 



T^-" = {e + P) u^u'' - Pg'"'. 



(3.2) 



m'^ = 7(1,v). 



(3.3) 



7 = (1 -t;2) 2 . 



(3.4) 
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In general situations, in the relativistic systems the energy density and pressure 
depend on the temperature and baryon number density. This requires that the 
baryon number conservation law, 

d,J^ = 0, (3.5) 

should be considered together with fl3.ip . However, in the case of the central rapidity 
region of ultra-relativistic heavy-ion collisions the dominating degrees of freedom are 
mesons (initially gluons), hence the net baryon number is very close to zero. In such 
a case the local value of baryon chemical potential is negligible, /ib ~ 0, and we can 
express the entropy density s, pressure P, and energy density e by temperature alone. 
This allows us to use the thermodynamic identities (12. ip and write the conservation 
laws (13.11) in the form. 

a^(TsM^M^) = 9^P. (3.6) 

After performing simple transformations shown explicitly in Appendix IB. 21 Eg. (13.61) 
leads to the two formulas 

u^d^iTu") = d'T, (3.7) 

d^isun = 0. (3.8) 

Eq. (13.71) is the acceleration equation and represents the relativistic analog of the 
Euler equation known from the classical hydrodynamics. Eq. (13. 8p states that the 
evolution is adiabatic (entropy is conserved). In the non-covariant notation the 
hydrodynamic equations (13. 7p and (13.81) are expressed by the following expressions 



^(T7v) + V(T7) = V X [V X (T7V)] , (3.9) 

^(S7)+V(s7v) = 0. (3.10) 

The four equations above can be written in the Cartesian coordinates in the equiv- 
alent form as 



2 , i9 In T dlnT 1 dv 

dt + dt ^l-v-'dt 


= 0, 


(3.11) 


2, / dlnT dlnT\ dv^ dvy 
' "H'^ dx dy J^'^^dt '^dt 


= 0, 


(3.12) 


o.dlnT dlnT dvz Vz dv 
dz dt + dt+l-v^'^dt 


= 0, 


(3.13) 


dins V dv dvx dvy dvz 
dt 1 — f ^ dt dx dy dz 


= 0, 


(3.14) 



where the total time derivative is defined by the equation 

d d 



^ , V - V. (3.15) 

dt dt ^ ^ 
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Figure 3.1: Decomposition of the flow velocity vector in the plane z = 0. In our approach 
we use the magnitude of the transverse flow vj_ = + v"^ and the angle a as two inde- 
pendent quantities, isted of Vx and Vy. The longitudinal component of the flow Vz = z/t, 
a consequence of boost-invariance. 



The details of the transformations leading from (13. 6p to (13. lip - (I3.14p can be found 
in Appendix IB. 31 Equations (13. lip - (13.140 do not form a closed system of equations 
since they contain five independent variables T, s, v^, Vy and Vz- An additional 
equation is needed to close them, i.e., the equation of state is required which in- 
troduces the relation between T and s. Alternatively, the equation of state may be 
included by the use of the temperature dependent sound velocity 

, dP sdT 

Cs(T) = ^ = ^^, (3.16) 

The temperature dependent sound velocity for strongly interacting matter that is 
used in our work was discussed thoroughly in Chapt. [2] and is plotted in Fig. 12. 1[ 



3.2 Implementation of boost-invariance 

The experimental data collected at RHIC by the BRAHMS Collaboration sug- 
gests that in the midrapidity region the particle yields do not vary much with rapid- 
ity. Thus, we can assume that number of particles per unit rapidity in the range of 
\y\ < 1 is essentially constant and the midrapidity region (central region) is boost- 
invariant. This symmetry demands that the longitudinal component of the velocity 
has the form of the Bjorken flow ^7\ . 



(3.17) 
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and the thermodynamic scalar variables like temperature or entropy density are 
functions of the longitudinal proper time r = ^t'^ — and the transverse coordi- 
nates X and y. In 'practice, these properties mean that we can solve the hydrodynamic 
equations for 2 = and by using the appropriate Lorentz transformations we obtain 
solutions for z 0. 

Adopting the procedure outlined above and restricting our considerations to the 
plane 2; = we observe that Eq. fl3.13p is automatically fulfilled and we are left with 
only three independent equations 

, n , 9 In T dlnT 1 dvi 

- ("»^ - "'^) + '"''li - ■"'IT = P-^^' 

dins v± dv± dvx dvy 1 ^ ^ 
dt 1 — v'j_ dt dx dy t ' 



where v± = + f ^ is the transverse velocity. The equations describing transverse 
evolution can be rewritten in the cylindrical coordinates which are convenient for 
further analysis [H] (details can be found in Appendix IB.4.11 and IB.S.ip 



2 c^lnT cilnT 1 dv^ 



(3.21) 



9 , / 9 In T v\ cos a 5 In T \ r, f da f 1 sin a \ 

dins 1 dv_[_ dv± 9 da cos a /da \ v± 
v±^r- + 2— — -vlsma— + — ( — + 1 ) + — = 0. 



(3.22) 



dt 1 — V , dt dt dr r \d(j) J t 

(3.23) 



Here r is the distance from the beam axis, r = a/x^ + y^, and is the azimuthal 
angle, = tan~^(?//x). These two coordinates parameterize the plane z = 0. The 
angle a is the function describing direction of the flow, a = tan~^{vy/vx) — 0, see 
Fig. 13.11 The differential operator d/dt used in Eqs. (13.211) - (13.231) is defined by the 
formula 

d d d v±sina d a oa\ 

dt dt ""^^ ^ dr r dd) 



3.3 Characteristic form of hydrodynamic 
equations 

In order to obtain the form of the hydrodynamic equations which is convenient for 
numerical studies we introduce new independent variables. In this respect we follow 
the method originally proposed by Baym et al. in Ref. [13] and the new form of 
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the equations will be called the characteristic form. We introduce the potential $ 
defined as 

d<!> = -dlnT = Csdlns, (3.25) 

Cs 

and the transverse rapidity r]± = ta,ioh~^ v±. The new form of the hydrodynamic 
equations is expressed by the dimensionless auxiliary functions and A_ defined 
by the equations 

A± = ^±ri±. (3.26) 

The sum and difference of Eq. f l3.2ip and fl3.23p together with Eq. fl3.22p take the 
form 



dAj 
~dt 



da 
'dt 



1 ± CsV± 
Cs 

1 ± CsV± 



dA+ sin a dA-t 



± bill a u^± 
COS a — 1 ;-— 



dr 



da v± cos a f da 



+ v± cos a — — h 
dr 



v± sm a — 

dr 

da v± sin a f da 



0, (3.27) 



Csfl 



7^ + 1 



9$ COS a 5$ 



sm a 



dr 



0, 



(3.28) 



where the transverse velocity and the $ potential are expressed through the auxiliary 
functions A± as follows 



v± = tanh 



A, 



A_ 



A.+A. 



(3.29) 



Temperature and all other temperature dependent variables, e.g., the sound velocity, 
can be also calculated from the functions A±, 



T = T<i,(<l>) =T$ 



A, 



Cs(T) 



2 

A^ + A_ 



(3.30) 
(3.31) 



Here we have introduced the subscripts to make clear what kind of the argument is 
expected for a given function. For example, the temperature may be considered as a 
function of entropy density or $. In those two cases one should use the functions Ts 
or T$, respectively. If the equation of state is known all such functions can be easily 
calculated and Eqs. (13.271) and f l3.28p may be used to determine three unknown 
functions A^, A_, and a. 



3.4 Boundary conditions 

Having in mind the heavy-ion collisions at RHIC and at LHC we consider the col- 
lisions of identical nuclei, Au+Au or Pb+Pb, which collide moving initially along 
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Figure 3.2: Non-central collision of two identical nuclei at the impact vector b, viewed 
in the transverse plane. The nucleus A is located at (xi,yi) = (—6/2,0) and the nucleus 
B at (x2,y2) = (6/2,0). The overlapping region has the shape of an almond elongated in 
the direction of the y-axis. 



the z-axis. The positions of the centers of nuclei depend on the impact parame- 
ter b and for non-central collision they may be located in the transverse plane at 
xi = {xi,yi) = (—6/2,0) and at X2 = {x2,y2) = {b/2,0), see Fig. 13.21 The distri- 
bution of matter created just after the collision is not cylindrically symmetric and 
has an ellipsoidal shape (one commonly speaks of an almond shape). With the use 
of the coordinate system explained above the transverse velocity must vanish at the 
origin of the system namely 

t;x(t,r = 0,0) = 0. (3.32) 

Also the gradients with respect to the distance r of temperature and entropy density 
must converge to zero at r = 

or or 

The boundary conditions (13.321) and (13.331) can be naturally fulfilled by the following 

Ansatz 

A+{t,r,(f)) = A(t,r,0), 

A_(t,r,0) = A{t,-r,(f)), r>0 (3.34) 
a{t,—r,(f)) = a{t,r,(l)). 

The domain of the transverse distance r has been extended to the negative values 
of r. The two functions A^ and A_ are replaced in this way by a single function A. 
At the same time the function a has been generalized to have negative arguments 
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by the condition that it is a symmetric function of r, see Fig. 13.31 With the help 
of the definitions (13.341) . Eqs. fl3.27p may be reduced to a single equation for the 
function r, (/>), 



dA ^ v± + Cs 

dt 1 + Cgf _L 



dA sin a dA 
cos a — H — 

or r od) 



1 + Csf ± 



v±_ sm a 



da f _L cos a f da 



dr 



0. (3.35) 



The transverse velocity and $ potential from Eq. (13.291) are now expressed as 

'Ait,r, 0)-A(t, -rj 



v±{t,r, (f)) 
$(t,r, 0) 



tanh 



A{t, r, 0) + A{t, -r, 0) 



(3.36) 
(3.37) 



In addition, the use of cylindrical coordinates implies the periodicity of all functions 
in angle thus creating another set of periodic boundary conditions 



A(t,r,0 = O) = A(t, r, = 27r) 
Q;(t, r, = 0) = a(t,r, = 27r) 



(3.38) 




-20 -10 10 20 

r [fm] 



Figure 3.3: Construction of the functions A±(t., r, tp) in terms of a single function A{t, r, 0), 
see Eq. (I3.34p . The function a(t, r, 0) for negative values of r is obtained from the 
symmetry condition a{t, —r, 0) = a{t, r, 0). 
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Chapter 4 
Initial conditions 



Along with the formulation of the hydro dynamic equations that govern the evolu- 
tion of matter we must also specify a set of initial conditions that are required to 
unambiguously solve such equations - the hydrodynamic equations are first-order 
partial differential equations. The physical systems studied here are formed in the 
collisions of two identical gold nuclei (RHIC experiment) or two identical lead nuclei 
(the future heavy- ion program at LHC). In this Chapter we discuss in more detail 
the initial conditions used by us to analyze the collisions at RHIC and LHC. At 
first we discuss the most common form of the initial conditions that is based on the 
Glauber model. Next, we present the modified initial conditions where the initial 
energy density has a form of the two-dimensional Gaussian in the transverse plane. 
Finally, we present our method of including the parton free-streaming as a process 
which precedes the hydrodynamic evolution. In particular, we show how to match 
the free-streaming stage with the hydrodynamic stage. 

4.1 Standard initial conditions 

In the calculations presented in Chapter [6] and [7] we assume that the initial entropy 
density of the particles produced at the transverse position point is proportional 
to a source profile obtained from the Glauber approach. This is done in the similar 
way to other hydrodynamic calculations. Specifically, for the particle source profile 
we use a mixed model ^91 [50] . with a linear combination of the wounded-nucleon 
density pwN and the density of binary collisions pec namely 

1 — K 

s(x_l) oc p(xi) = — ^Pwn(x_l) + kpbc(x±). (4.1) 

The case k = corresponds to the standard wounded-nucleon model [HI], while 
K = 1 would include the binary collisions only. The PHOBOS analysis [50] of 
the particle multiphcities yields k = 0.12 at ^snn = 17 GeV and k = 0.14 at 
y^SNN = 200 GeV. Since the density profile from the binary collisions is steeper 
than from the wounded nucleons, increased values of k yield steeper density profiles, 
which in turn result in steeper temperature profiles. 
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In our hydro dynamic code, the initial conditions are specified for the temperature 
profile which takes the form 



Tin, x_l) = Ts 



P(0) 



(4.2) 



where Ts{s) is the inverse function to the function s{T), and Si is the initial en- 
tropy at the center of the system. The initial central temperature Ti equals Ts{si). 
Throughout this work, the initial time for the start of the hydrodynamic evolution 
is denoted by Xj. 

The wounded-nucleon and the binary-collisions densities in Eq. 04.21) are obtained 
from the optical limit of the Glauber model, which is a very good approximation for 
not too peripheral collisions [52]. The standard formulas are [51] 



PWN IX_lJ 



and 



Ta 



+ Ta 



PBC (XxJ 




■Ta 



0"in Jyl I — 



A 



^±]Ta{ -- 



-I A' 



-I A' 



(4.3) 



(4.4) 



In Eqs. (14.31) and (14.41) b is the impact vector, o"in is the nucleon-nucleon total 
inelastic cross section, and T4 (x, y) is the nucleus thickness function 



TAix,y) = / dzp{x,y,z) . 



(4.5) 



For RHIC energies we use the value = 42 mb, while for LHC we take (Tin = 63 mb. 
The function p{r) in Eq. (14.51) is the nuclear density profile given by the Woods- 
Saxon function with the conventional choice of the parameters: 

po = 0.17fm-^ 

ro = (1.12A^/3-0.86v4-^/^)fm, 

a = 0.54 fm. (4.6) 

The values of the atomic mass A are: 197 for RHIC (gold nuclei) and 208 for LHC 
(lead nuclei). The value of the impact parameter in Eqs. (14. 3 p and (14.41) depends 
on the considered centrality class. 

Besides the initial temperature profile (14. 2 p we also specify the initial transverse 
fiow profile. 



v{Ti,r,, 



Hor 



0. 



(4.7) 
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The results presented below are obtained with Hq = 0.001 fm^^. The very small 
value of the parameter Hq means that the system is practically at rest at the mo- 
ment when the hydrodynamic evolution starts. However, nonzero Hq improves the 
stability of the numerical method. 

In our numerical calculations we use the adaptive method of lines in the way 
as it is implemented in the MATHEMATICA package. The and r directions are 
discretized, and the integration in time is treated as solving of the system of ordinary 
differential equations. Typically, we use grids with Ar = 0.25 fm in the r direction 
and A(j) = 6 degrees in the (p direction. This method fails in the case where the 
shock waves are formed, however with our regular equation of state and the regular 
initial conditions such shocks are not present. 

We stress that the shape of the initial condition (14. 2 p is important, as it deter- 
mines the development of the radial and elliptic flow, thus affecting such observables 
as the pT-spectra, V2, and the femtoscopic features. On qualitative grounds, sharper 
profiles lead to more rapid expansion. Several effects should be considered here. 
Firstly, as discussed in Ref. hydrodynamics may start a bit later, when the pro- 
file is less eccentric than originally due to initial free-streaming of partons in the 
pre-hydro phase. On the other hand, statistical fluctuations in the distribution of 
the Glauber sources (wounded nucleons, binary collisions) |53l IMl EH [561 EH EH ES] 
lead to a significant enhancement of the eccentricity, especially at low values of the 
impact parameter. Thus the initial eccentricity may in fact be smaller or larger than 
what follows from the application of the Glauber model. This contributes to the 
systematic model uncertainty at the level of about 10-20%. This uncertainty could 
only be reduced by the employment of a realistic model of the pre-hydro dynamic 
evolution. With this uncertainty in place, one should not expect or demand a better 
agreement with the physical observables than at the corresponding level of 10-20%. 

The results obtained with the introduced here standard initial conditions 
are presented and discussed in Chapters [H] and 

4.2 Gaussian initial conditions 

The situation described in the previous Section corresponds to the typical case 
where the hydrodynamic evolution is initiated from a source profile generated by 
the Glauber model, with the initial central temperature or energy serving as a 
free parameter. The initial density and possibly the initial flow profiles may be also 
provided by the early partonic dynamics, for instance by the Color Glass Condensate 
(CGC) [nUl [M] . In practice, however, the theory of the partonic stage carries some 
uncertainty in its parameters, which influences our knowledge of the details of early 
dynamics. Having in mind such uncertainties we decided to try another class of the 
initial conditions, i.e., in our calculations we also include the case where the initial 
energy profiles in the transverse plane have the Gaussian shape. 

The modified Gaussian parameterization of the initial energy density profile at 
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the initial proper time t; has the following form 



£(x^j = Ei exp 



2^ ~ 262 



(4.8) 



The energy density profile (14.81) determines the initial temperature profile that is 
used in the hydrodynamic code, 



7'(ri,x_L) = 



Ei exp 



?/2 



2a2 262 



(4.9) 



The initial central temperature T(ri,0), which may depend on the centrality, is de- 
noted by Ti and is a free parameter of our approach. 



The results obtained with such initial conditions are presented and dis- 
cussed in Sect. 18. IL 



4.3 Free streaming 

The very early start of hydrodynamics means that the matter equilibrates very fast. 
Such a short thermalization time is difficult to explain on the microscopical grounds 
and inspires hot discussion about the nature of matter produced at RHIC. To avoid 
the problem with the early thermalization we consider also the scenario where the 
hydrodynamic evolution is preceded by the free-streaming stage. In this version 
of our calculations we assume that partons behave as free particles in the proper 
time interval 0.25 < r < 1 fm. At = 1 fm the sudden equilibration takes place 
which is described with the help of the Landau matching conditions. The global 
picture is as follows: early phase (CGC) generating partons at time tq = 0.25 fm, 
partonic free streaming until Ti = 1 fm, hydrodynamic evolution until the freeze-out 
at temperature Tf , free streaming of hadrons and decay of resonances. We note that 
similar ideas have been described by Sinyukov et al. in Refs. [021 E3]- 

In the remaining part of this Section we describe in more detail the free-streaming 
stage. We assume that massless partons are formed at the initial proper time 
To = \/to — and move along straight lines at the speed of light until the proper 
time when free streaming ends, r = \fW^^^z^ . We introduce the space-time rapidi- 
ties 

,„ = ilog^ (4.10) 

and 

')ll = 5log^. (4.11) 

Elementary kinematics ^62j links the positions of a parton on the initial and final 
hypersurfaces and its four-momentum 

= (pt cosh y, pt cos 0p, px sin 0p, px sinh y), (4-12) 
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where y and px are the parton's rapidity and transverse momentum: 

Tsinh(?7|| -y) =To sinh(?7o - y), 

X = xo + A cos (pp, y = yo + A sin (^p, 

A = T cosh(?/ — rjw) — \ Tq + sinh^(?/ — ?7||). 



(4.13) 



Thus the phase-space density of partons at the proper times tq and r are related 



dy(PpTdrj\\dxdy 



drjodxadyo 



d^N 



S { Vo ~ y — arcsinh 



dyd^PTdrjodx^dyQ 

T 



X 



^0 



sinh(?7|| — y) 



(4.14) 



A sin I 



S{x-Xo- A cos (l)p)S{y-yo ^^...^p, 

In our approach we restrict ourselves to the boost-invariant systems, hence it is rea- 
sonable to assume the following factorized form of the initial distribution of partons, 

d^N 

— — — , , , = n{xo,yo)F{y -r]o,pT), (4.15) 
dyd'^PTdriQdxodyo 

where n is the transverse density of partons obtained again from the GLISSANDO 
model, namely 



n{xo,yo) = noexp 



Xn 



2a2 262 



(4.16) 



When the rapidity emission profile F is focused near y = rjQ, for instance if we have 
F ~ exp[(?/ — rjo)'^/ (2 6y'^)], with 6y ~ 1, and if r ^ tq, then the kinematic condition 
(14.131) effectively transforms it into 



arcsinh 



exp 



^ sm{y - m 



2 5y^ 



(4.17) 



Then Eq. Km yields 

d^N 



n{x- At cos (f)p,y- At sin (t>p)6{y - r]\\)f{pT) 



(4.18) 



dyd'^PTdr]\\dxdy 

where At = t — tq and fipr) is the transverse momentum distribution. The energy- 
momentum tensor at the proper time r, rapidity r^y , and transverse position (x, y) 
is given by the formula 

d^N 

dyd^PT-f-p^—r—f-rP'^p" (4.19) 



rpflU 



2lT 



dyd^prdriwdxdy 



,n{x — Ar cos0p, y 



At sin (f)p) x 



cosh ?7|j cos (pp cosh ?7|j sin 0p cosh r^y sinh r^y \ 



27mo 
/ cosh^ ?7|| 
cosh r^ii cos ( 
cosh sin ( 

y cosh sinh 7]^\ cos (j)p sinh ri\\ sin (f)p sinh ; 



cos 



cos (pp sm (pp 



cos (pp sm (pp 



sm 



COS( 



, sinh ?7|| 
sin 0p sinh rj\\ 
sinh^ T^ll 
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Figure 4.1: Sections of the energy-density profile e (Gaussian-like curves) normalized to 

unity at the origin, and of the velocity profile v± = -^Jv"^ + Vy (curves starting at the 

origin), cut along the x axis (solid lines) and y-axis (dashed lines). The initial profile is 
from Eq. (j4.16p for centrality 20-40% at r = tq. The e profiles are for r = tq = 0.25, 1, 
and 2 fm, and the velocity profiles are for r = 0.25, 1 and 2 fm, all from bottom to top. 
We note that the flow is azimuthally asymmetric and stronger along the x axis. 



where the factor eo/(2v™o) is a constant from the px integration Due to boost 
invariance the further calculation may be carried for simplicity of notation at r)\\ = 0. 
Next, we assume that at the proper time r the system equilibrates rapidly. We thus 
use the Landau matching condition, 

T^'^(x, y)u,{x, y) = e{x, y)g^''u,{x, y), (4.20) 

to determine the position dependent four-velocity of the fluid, 

M'^(r, X, y) = — i=(l, v^, Vy, 0) (4.21) 

and its position dependent energy density e{T,x,y), which is then identified with 
the energy-density profile. The obtained in this way the flow and energy profiles are 
used as an input for the hydrodynamic calculations. The time r coincides with the 
starting time for hydrodynamics, t = T[, and we again use the formula 

T(ri,xx) =T,[£(ri,0,0)] (4.22) 

to define the initial temperature profile in the hydrodynamic code. The initial 
central temperature T^ = Ti,{ei) is again used as a free parameter (connected with 
the freedom of choosing eo). 

The results of solving Eq. fH:20|) with Tf"" from (Hl9D for Ar = 0.75 fm and the 
initial profile (14.161) are shown in Fig. 14. 1[ The presented results correspond to the 
centrality class 20-40%. The curves normalized to unity at the origin, r = 0, show 
the sections along the x and y axes of the energy-density profile at the proper times tq 

^Note that the structure of Eq. (|4.19p implies that for A = the initial energy density given 
by T*"^ and the initial parton density n have the same profiles in the transverse plane. 
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(no free streaming) and r. Obviously, at tq we find e{x,y) = eon{x,y)/no. We note 
that the profile spreads out as the time progresses. Importantly, this effect is faster 
along the shorter This is clearly indicated by the velocity profiles along the 

X and y axes (the curves starting from at r = 0). Thus the fiow generated by free 
streaming and sudden equilibration is azimuthally asymmetric. This asymmetry 
is also seen from the expansion in the parameters At, x,y, where straightforward 
algebra gives (for xAt « a? and yAr « b"^) 




(4.23) 



The results of the hydrodynamic calculations preceded by the free-streaming 
stage are presented and discussed in Section 18.21 
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Chapter 5 

Freeze-out prescription 



As the system expands, its density decreases and the mean free path of particles 
becomes larger and larger. When the average time between the collisions becomes 
compatible with the characteristic time of the system expansion, that is given ap- 
proximately by the inverse of the divergence of four- velocity, Texp ~ the 
hadrons decouple and may be regarded as free particles. The hydrodynamic stage of 
the evolution of matter is then changed into the free streaming of hadrons which are 
later registered by various detectors. The transition from strongly interacting hydro- 
dynamic system into the non-interacting system of hadrons is called the freeze-out 
process. In fact one can distinguish two types of freeze-outs: the chemical freeze-out 
(when the inelastic processes stop and the hadronic abundances are frozen) and the 
kinetic freeze-out (when all interactions stop and the momentum spectra are frozen). 
Of course, the kinetic freeze-out may be regarded as the true freeze-out, after which 
all interactions cease. The concept of the chemical freeze-out is however very use- 
ful since this is the stage when the hadronic abundances are fixed (by definition, 
the subsequent evolution is dominated by elastic collisions which do not change the 
hadronic yields). 

In our approach we use the simplification that the two freeze-outs discussed 
above happen at the same time. This approach is known as the single-freeze-out 
assumption and several earlier studies showed that it leads to the correct description 
of the data [30l [6ll [65l [661 EH EH] . Of course in the hydrodynamic approach the 
freeze-out happens locally, i.e., at different times for different space positions. In our 
studies we assume that freeze-out takes place at a given final temperature Tf . At the 
beginning of the evolution the hadrons are emitted usually from the edge (surface) 
of the system where the temperature is always low. At the end of the evolution, the 
hadrons are emitted mainly from the whole volume, since the temperature eventually 
drops down to Tf everywhere due to the hydrodynamic expansion. 

5.1 Cooper- Frye formula 

The collection of the spacetime points where the freeze-out process takes place is 
called the freeze-out hypersurface. This hypersurface is a three dimensional manifold 



43 



44 



CHAPTER 5. FREEZE-OUT PRESCRIPTION 




Figure 5.1: Parameterization of the freeze-out hyper surf ace. The part (a) represents the 
view in the r — r plane with the fixed values of and r/|| . The part (b) shows the view in 
the X — y plane with the fixed values of ( and r^y . 

in the four-dimensional Minkowski space. The approach which is used to calculate 
the spectra of the emitted hadrons is based on the Cooper-Frye formula [691 EO] 

-j^= [ dE,p^f{p-u). (5.1) 
dydy± J 

In Eq. (15.1 p dT,^ is the element of the freeze-out hypersurface and / is the equilibrium 
distribution function depending on the scalar product of the particle four-momentum 
and the fluid four- velocity u'^. In our approach the form of the freeze-out hy- 
persurface is delivered by the hydrodynamic code. The criterion T = Tf is used 
to find the appropriate space-time points, whose positions are parameterized in the 
following general way 

t = [Ti + d (0, C, sin C] cosh r/n , z= [n + d (0, C, V\\) sin C] sinh ^ , 

X = (0,C,77||) cosC COS0, y = d{(j),C,V\\) cosC sm(j). (5.2) 

Here is the standard azimuthal angle in the y — x plane, while the variable m is 
the space-time rapidity defined as 

ry|| = -ln^ = tanh- (-j . (5.3) 

The parameterization (15.21) yields the compact expressions for the proper time r and 
the transverse distance r. 



T-Ti = Vt"^ - z"^ - Ti = d {(f), Cm) sinC, 

r = + y2 = d (0,C,'7||) cosC- (5.4) 
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For fixed (f) and ?7|| the quantity d (0, (, T^y) is the distance between the hypersurface 
point with coordinates (0, C,?7||) and the spacetime point (r = Xj, x = 0, y = 0), see 
Fig. 15. 1[ The variable (, restricted to the range <(< it /2, is an angle in the 
T — r space. We have introduced the angle ( because in most cases the freeze-out 
curves in the t — r plane may be treated as functions of this parameter. The use 
of the transverse distance r is inconvenient, since very often two freeze-out points 
correspond to one value of r. 

If the parameterization of the freeze-out hypersurface in terms of the variables 0, 
C, and 7711 is given, then the element of the freeze-out hypersurface may be obtained 
from the formula known from the differential geometry 



dT^a 



-a/375 



d4> 

dx-t 
d<t> 



dx^ dxl^ 



dx-f 
d( 
dx^ 



dx]_ 

dm 
a»7|| 



dcj) dC (i?7|| , 



(5.5) 



where eajj-ys is the totally antisymmetric tensor with eoi23 = +1- The direct differ- 
entiation gives: 



dT,^ = d cos ( 



(ti + d sin Q (^d sin ^ — |^ cos cosh r^y + ^ sinh 1 



dE^ = -d (Ti + rf sinC) 
dS^ = -d (t; + d sin C) 



d cos C + sin C j cos C cos + |^ sin i 



d cos C + sin C j cos C sin + || cos 1 



(5.6) 



dS^ = d cos C 



(xi + d sin () (^d sin ^ — |^ cos sinh r^n + ^ cosh 1 



With Eqs. fl5.6p and the standard definition of the four-momentum expressed in 
terms of the rapidity and transverse momentum. 



= {rriT cosh y, px cos 0p, px sin (pp, uit sinh y) 



(5.7) 



where my = a/ rn^ + p\ is the transverse mass, we find the explicit form of the 
Cooper-Frye integration measure 



dT,^p'^ = d{Ti + dsin() 
dd 



d cos ( {niT sin ^ cosh (?7|| — y) + Pt cos C cos (0 — 0p)) 



-|- — cos C {—nriT cos C cosh (?7|| — y) + Pr sin C cos (0 — 0p)) 



dd dd > . , / \ 

+ -q^Pt sm (0 - 0p) + ^"^T- cot C smh (r/y - y ) 



dri\\d(f)d(. 



(5. 
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This equation, when used in the Cooper- Frye formula [HS], leads to the six- dimen- 
sional particle density at freeze-out 

dN 

dyd(j)pPTdpTdrj\\ d(j)dC. 



In the transition from (15.81) to (15.91) we used the assumption of boost-invariance and 
removed the term dd/dri\\ - for boost-invariant systems the function d depends only 
on and C,. The last line in (15. 9p contains the equilibrium statistical distribution 
function, with —1 for bosons and +1 for fermions. The argument is the Lorentz- 
invariant product p^u^. The parameter (3 = 1/T is the inverse temperature and g 
denotes the spin degeneracy factor. The four-velocity field has been expressed in 
terms of the transverse flow v±{(j), () and the dynamical angle a(0, () which depend 
on the space-time positions on the freeze-out hypersurface. 

Although in the hydrodynamic evolution we neglected the effects of the chemical 
potential fi, it is introduced in the generation of the particles at freeze-out. Thus, 
the equilibrium distributions include the term 



with B, S, and /a denoting the baryon number, strangeness, and isospin of the 
particle. The baryonic, strange, and isospin chemical potentials assume the values: 

= 28.5 MeV, /is = 9 MeV, /ijg = -0.9 MeV at the highest RHIC energies 
and fiB = 0.8 MeV and /xg = fii^ = MeV at the LHC [33]. 



5.2 Calculation of observables 

In order to calculate the physical observables we first run our hydrodynamic code and 
extract the functions describing the freeze-out hypersurface and flow: d = d{(f),()j 
v± = C)) and a = a{(j),C)- In the next step these functions are used as 

input for the THERMINATDR|27j code, which generates, according to the formula (15. 9p . 
the distributions of the primordial particles. The primordial particles include the 
stable hadrons as well as all hadronic resonances. The resonances decay through 




(5.9) 



(5.10) 
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strong, electromagnetic, or weak interactions at the random proper time controlled 
by the particle life-time, and at the location following from the kinematics. As the 
final result we obtain the stable particle distributions including the feeding from 
the resonance decays. Since the memory on all decays is kept, one may apply 
the experimental cuts or the weak-decay feeding policy, which facilitates the more 
accurate comparison of the model to the data. For example, one may extract easily 
the model proton spectrum which does not include the feeding from the A decays. 
The momentum distribution functions are obtained by integrating over 0, ri\\, and 

c, 

oo 2tt 7r/2 

dN r , r , r dN 



dyd^PT J J ^ dyd(l)pPTdpTdri\\d(l)d( ^ ^ 

-oo 

For cylindrically asymmetric collisions and midrapidity, y = 0, the transverse- 
momentum spectrum has the following expansion in the azimuthal angle of the 
emitted particles 

dN dN , / N 

:i + 2t;2(pT)cos(20p) + ...). (5.12) 



dyd?j)T dy 27rpT dpT 

Equation f l5.12p defines the elliptic flow coefficient V2- 

The well-known problem of the freeze-out description in the Cooper- Frye formu- 
lation is that the hypersurface from hydrodynamics typically contains a non-causal 
piece, where particles are emitted back to the hydrodynamic region [7T]. In our ap- 
proach we follow the usual strategy of including only that part of the hypersurface 
where dY^^p^ > 0. To estimate the effect from the non-causal part we compute the 
ratio of particles flowing backwards {i.e. where dT^^p^ < 0) to all particles. For the 
cases studied in the following sections we find this ratio to be a fraction of a percent. 
Such very small values show that the known conceptual problem is not of practical 
importance for our study. 

The correlation function for identical pions is obtained with the two-particle 
Monte-Carlo method discussed in detail in Ref. [721 173] . In this approach the eval- 
uation of the correlation function is reduced to the calculation of the following 
expression 

EE^A(q-p. + P,) 5A(k-i[pi + p,]) |vl/(k*,r*)|2 

= ^ '''VvA ( ^ w M ^ n ' (^-^3) 

L L <^A (q - Pi + Pi) ()a (k - 2 [p* + Pj\) 
where 5^ denotes the box function 

if bxl < f,bj < f < f ^g^^^^ 

otherwise. 




0, 



In the numerator of Eq. 05.131) we include the sum of the squares of modules of the 
wave function calculated for all pion pairs with the relative momentum q (we use 
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the bin resolution A = 5 MeV) and the pair average momentum k. For non-central 
collisions we only provide azimuthally integrated HBT radii. For each pair the wave 
function \E'(k*, r*), including the Coulomb interaction, is calculated in the rest frame 
of the pair; k* and r* denote the relative momentum and the relative distance in the 
pair rest frame, respectively. In the denominator of Eq. f l5.13p we put the number of 
pairs with the relative momentum q and the average momentum k. The correlation 
function fl5.13p is then expressed with the help of the Bertsch-Pratt coordinates 
kr, Qout, Q'side, Qiong and approximated by the Bowler-Sinyukov formula [HI [75] 

C(q, k) = (1 - A) + A K^oniiq-mv) X 

X [l + exp ( — -Rout Q'out " -^side Iside ~ -^long ftong)] ' (5.15) 

where -ft'couK^inv) with ginv = 2fc* is the squared Coulomb wave function integrated 
over a static Gaussian source. We use, following the STAR procedure [TB], the 
static Gaussian source characterized by the widths of 5 fm in all three directions. 
Four kr bins, (0.15 - 0.25), (0.25 - 0.35), (0.35 - 0.45), and (0.45 - 0.60) GeV, 
are considered. The 3-dimensional correlation function with the exact treatment of 
the Coulomb interaction is then fitted with this approximate formula and the HBT 
radii -Rout, -Rside, and -Riong are obtained. They can be compared directly to the 
experimental radii. 



Chapter 6 

Soft-hadronic observables at RHIC 



In this Section we apply our model to describe soft hadron production in the rel- 
ativistic heavy-ion collisions at the highest RHIC energy, y/s^N = 200 GeV. We 
analyze the data dehvered by the PHENIX [TTJ [7H] and STAR Collaborations flE\. 
We consider the most central events, defined by the centrality class c = - 5%, 
where we use the impact parameter b = 2.26 fm. We also use the data from the 
centrality classes c = 20 - 30% (transverse-momentum spectra and the HBT radii) 
and c = 20 - 40% (elliptic fiow), for which we take (for simplicity) the same value 
of the impact parameter, b = 7.16 fm. 

The link between the centrality class and the impact parameter involves the 
total Au+Au inelastic cross section. To a very good accuracy for not-too-peripheral 
collisions we may use the formula 



(c) ^ 3- (6-1) 

in 

The total inelastic Au+Au cross section crj^""^" is not directly measured. Typically, 
it is obtained from Glauber simulations and carries model uncertainties [52] • With 
total inelastic pp cross section (Tin = 42 mb and for gold nuclei used at RHIC, the 
GLISSANDO Monte-Carlo Glauber simulation [57] gives cr^^"^" ~ 6.4 b, so 6 = 2.26 
fm for (c) = 2.5% and b = 7.16 fm for (c) = 25%. 



6.1 Central collisions 

First, we consider the centrality class - 5% with the corresponding value of the im- 
pact parameter b = 2.26 fm. The in-plane and out-of-plane freeze-out curves are de- 
fined as the intersections of the freeze-out hypersurface with the planes y = ((^ = 0") 
and X = 0{(f) = 90°). They are obtained with the initial starting time for hydrody- 
namics Ti = 1 fm, the initial central temperature Ti = 320 MeV, the final (freeze-out) 
temperature Tf = 150 MeV, and are shown in Fig. 16. 1[ The two freeze-out curves 
practically overlap, indicating that the expansion of the system is almost azimuthally 
symmetric in the transverse plane. This feature is certainly expected for the almost 
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Figure 6.1: In-plane and out-of-plane freeze-out curves, i.e., tlie intersections of the freeze- 
out hypersurface with the planes y = 0{(p = 0°) and x = 0((j) = 90°), obtained for central 
RHIC collisions; c = - 5%, 6 = 2.26 fm, n = 1.0 fm, Ti = 320 MeV, and Tf = 150 MeV. 
The two curves overlap, indicating that the system at freeze-out is almost azimuthally 
symmetric in the transverse plane. 
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Figure 6.2: Transverse-momentum spectra of tt"'", K~^, and protons. The PHENIX experi- 
mental results [77] for Au+Au collisions at ^snn = 200 GeV and the centrality class 0-5% 
(points) are compared to the model calculations (solid lines) with the same parameters as 
in Fig. 16. 1[ 



central collisions, where the impact parameter is very small. We note that the shape 
of the isotherms is consistent with the result presented in Fig. 2 of Ref. [HO], where 
the same decoupling temperature of 150 MeV was considered. 

In Fig. 16.21 we present our results for the hadron transverse-momentum spectra 
with the same values of the parameters. The dots show the PHENIX data [77j for 
positive pions, positive kaons, and protons, while the solid lines show the results 
obtained from our hydrodynamic code linked to THERMINAT0R|27j . Our model de- 
scribes the data well up to the transverse-momentum values of about 1.5 GeV. For 
larger values of pt the model underpredicts the data. This effect may be explained 
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Figure 6.3: The pionic HBT radii plotted as functions of the average transverse momen- 
tum of the pair compared to the STAR data [76] at the centrality 0-5%. The calculation 
uses the two-particle method and includes the Coulomb effects. The values of the model 
parameters are the same as in Fig. 16.11 

by the presence of the semi-hard processes, not included in our approach. 

In our approach, the values of Ti and b control the overall normalization. On 
the other hand, the value of the freeze-out temperature Tf determines mainly the 
relative normalization and the slopes of the spectra. We note that the correct slope 
for the pions and kaons is recovered at a relatively high value of Tf . This is possible, 
since the spectra of the observed hadrons contain the contributions from all well 
established hadronic resonances. This single- freeze-out picture ^28l [Ml [65], where 
the same temperature is used to describe the ratios of hadronic abundances and 
the spectra, was first tested for hadronic spectra in Ref . [6l] , see also [HH [66] . We 
note that the value Tf = 150 MeV agrees with the recent results obtained in the 
framework of the single-freeze-out model in Ref. [H7]. In the calculation of the 
proton spectra, in order to be consistent with the PHENIX experimental procedure, 
we removed the protons coming from the A decays. 

In Fig. 16.31 the model results and the STAR data [7^ for the HBT radii are 
presented. Again, for the same values of the parameters, a quite reasonable agree- 
ment is found. Discrepancies at the level of 10-15% are observed in the behavior 
of the -Rside, which is too small, and -Riong, which is too large, probably due to the 
assumption of strict boost invariance. The ratio i?out /-Rside — 1-25 is larger than 
one, which is a typical discrepancy of hydrodynamic studies. Nevertheless, when 
compared to other hydrodynamic calculations, our ratio i?out/-Rside is significantly 
closer to the experimental value. The ratio is almost constant as a function of /c^, 
contrary to the decrease observed in the data. 
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Figure 6.4: The freeze-out curves for c = 20 - 30%, Ti = 1.0 fm, T, = 320 MeV, 
Tf = 150 MeV, and b = 7.16 fm. The sohd Hne describes the in-plane profile, while the 
dashed line describes the out-of-plane profile. 
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Figure 6.5: Transverse-momentum spectra of vr^, , and protons. The PHENIX ex- 
perimental results [77| for Au+Au collisions at ^/snn = 200 GeV and the centrality class 
20-30% (points) are compared to the model calculations (solid lines). The values of the 
model parameters are the same as in Fig. 16.41 

6.2 Non-central collisions 

Next, we consider the centrality classes 20 - 30% and 20 - 40%. The data are com- 
pared with the model results obtained with the impact parameter b = 7.16 fm. The 
values of the initial central temperature and the final temperature are the same as 
in the case of the central collisions, i.e., T; = 320 MeV and Tf = 150 MeV. The 
two freeze-out curves are shown in Fig. 16. 4[ We observe that the out-of-plane pro- 
file is wider than the in-plane profile. This difference indicates that the system is 
elongated along the y axis at the moment of freeze-out. This feature is in qualita- 
tive agreement with the HBT measurements of the azimuthal dependence of -Rside- 
The comparison of the experimental and model transverse-momentum spectra is 
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Figure 6.6: The elliptic flow coefficient V2- The values measured by PHENIX [78] at 
^/sNN = 200 GeV and the centrality class 20-40% are indicated by the upper (pions 
+ kaons) and lower (protons) points, with the horizontal bars indicating the px bin. 
The corresponding model calculations are indicated by the solid lines, with the bands 
displaying the statistical error of the Monte-Carlo method. The parameters are the same 
as in Fig. 16. 4i 
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Figure 6.7: The pionic HBT radii plotted as functions of the average transverse momen- 
tum of the pair compared to the STAR data [76] at centrality 20-30%. As in the case of the 
central collisions, the calculation uses the two-particle method and includes the Coulomb 
effects. The values of the model parameters are the same as in Fig. 16.41 

presented in Fig. 16. 5[ In this case we also find a reasonable agreement between the 
data and the model results. 

In Fig. 16.61 the model results and the PHENIX data [7S] on the elliptic flow 
coefficient V2{pt) are compared. We observe that the V2 of pions+kaons agrees with 
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the data. Taking into account the uncertainly in the initial eccentricity, discussed in 
Sect. 14.11 which is at the level of 10-20%, the obtained agreement is reasonable. On 
the other hand, the model prediction for proton V2 is too large. The discrepancy is 
probably caused by the final-state elastic interactions, not included in our approach 
[82]. 

The HBT radii for non-central RHIC collisions are shown in Fig. 16. 7[ Similarly 
to the central collisions, we observe that -Rside is slightly too small and -Riong is 
slightly too large. Still, the ratio i?out /-Rside is very close to the data. Comparing 
our values of -Rside with other hydrodynamic calculations we find that our values are 
larger. This effect is caused by the halo of decaying resonances which increases the 
system size by about 1 fm, see Ref. [72]. 



Chapter 7 

Predictions for LHC 



We expect that the increase of the initial beam energy from y/s^N = 200 GeV at 
RHIC to ^/sNN = 5.5 TeV at LHC resuhs essentially in a higher initial temperature 
Tj used as the input for the hydrodynamic calculation (in the midrapidity region 
studied here). Therefore, to make predictions for the collisions at the LHC energies 
we use a set of values for Ti which are higher than those used at RHIC, namely, 
Ti = 400, 450, and 500 MeV. Of course, it is not obvious which value of will be 
realized at LHC. Estimates based on extrapolations [83j suggest an increase of the 
multiplicity by a factor of 2 compared to the highest RHIC energies, which would 
favor Ti around 400 MeV. However, to investigate a broader range of possibilities, 
we take into account much higher temperatures as well. In the case of LHC we 
also use a larger value of the nucleon-nucleon cross section in the definition of the 
initial conditions (cXin = 63 mb, k, = 0.2, A = 208) and different values of the 
chemical potentials in the hadronization made by THERMINATOR (/iB = 0.8 MeV and 
fJ's = A^is = 0). The changes of ain and the nuclear profile imply that for LHC we 
use b = 2.4 fm for c = 2.5% and b = 7.6 fm for c = 25%. 

In Tables 17.11 and 17.21 we list our results for the following quantities: the total 
Tr"*" multiplicity dN/dy, the inverse-slope parameter for pions A, the pion+kaon el- 
liptic flow f 2 at = 1 GeV, and the three HBT radii calculated at the pion pair 
momentum fc^ = 300 MeV. The inverse-slope parameter given in Tables mi and \T2\ 
is obtained from the formula 



The tables show several expected qualitative features. Obviously, as the initial 
temperature increases, the multiplicity grows. This is due a larger initial entropy, 
which causes a larger size of the freeze-out hypersurface. We find that the following 
approximate parameterizations work very well for the multiplicity of vr^ at LHC for 
the two considered centrality cases: 




(7.1) 



dN 
dy 



12000 (Ti/GeV) 



3.4 



(6 = 2.4 fm) 



(7.2) 



dN 
dy 



6000 (Ti/GeV) 



3.4 



{b = 7.6 fm). 
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The power 3.4 works remarkably well. This behavior reflects the dependence of the 
initial entropy on T as shown in top left panel of Fig. 12.21 of Chapter [21 where for 
the relevant temperature range of 300 — 500 MeV we have the approximate scaling 
s/T^ ^ TO-l 

Similarly, for the slopes in the studied domain we have 

A(l GeV) = 0.639 T, + 0.033 GeV, (6 = 2.4 fm) 

(7.3) 

A(l GeV) = 0.661 + 0.033 GeV, (6 = 7.6 fm). 
The HBT radii increase rather moderately with T[, as can be seen from the tables. 



Ti [MeV] 


dN 
dy 


A [MeV] 


-Rsido [fm] 


-Rout [fm] 


-Rlong [fm] 


320 


274 


237 


4.2 


5.4 


6.5 


400 


543 


288 


5.2 


5.9 


7.7 


450 


802 


320 


5.6 


6.1 


8.2 


500 


1136 


352 


6.0 


6.2 


8.8 



Table 7.1: Central collisions at RHIC corresponding to the results from Chapter 16.11 (the 
second row) and LHC (the three lower rows): A set of our results obtained for four different 
values of the initial temperature: Ti = 320, 400, 450, and 500 MeV. The columns contain 
the following information: dN/dy - the total pion multiplicity (positive pions only), A - 
the inverse-slope parameter for positive pions at pT = 1 GeV, -Rside; -Rout) -Riong ~ the three 
HBT radii calculated at the average momentum kx = 300 MeV. 



Ti [MeV] 


dN 
dy 


A [MeV] 




-Rsidc [fm] 


-Rout [fm] 


-Rlong [fm] 


320 


152 


244 


0.14 


3.6 


4.5 


5.6 


400 


272 


299 


0.16 


4.2 


4.5 


6.2 


450 


401 


331 


0.16 


4.6 


4.6 


6.7 


500 


569 


363 


0.15 


5.2 


4.9 


7.4 



Table 7.2: Non-central collisions at RHIC corresponding to the results from Chapter 
16.21 (the second row) and LHC (the three lower rows): The same quantities shown as in 
Table [7T] with the additional information on the pion elliptic flow V2 at = 1 GeV. 



7.1 Central collisions 

In the following Sections we show more details, discussing the spectra, V2, and the 
HBT radii obtained for the initial temperature T^ = 400, 450 and 500 MeV. 
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Figure 7.1: The freeze-out curves for central LHC collisions with c = — 5%, b = 2.4 fm, 
Ti = 1.0 fm for three initial temperatures Tj = 400,450 and 500 MeV, and Tf = 150 MeV. 
Similarly to the central RHIC collisions, the two curves overlap indicating that the system 
at freeze-out is symmetric in the transverse plane. 




Figure 7.2: The model results for the transverse-momentum spectra of vr"'", K~^, and 
protons. The values of the model parameters are the same as in Fig. 17.11 



In Fig. 17.11 we show the freeze-out curves obtained from our hydrodynamic code 
with Ti = 400, 450 and 500 MeV. Comparing to the corresponding central RHIC 
collisions with Ti = 320 MeV from Fig. 16.11 we observe that the difference in the 
initial temperature results in the longer time of the hydrodynamic expansion and a 
larger transverse size (both increase by about 3 fm). On the other hand, similarly to 
the RHIC results, we find that the two freeze-out profiles overlap, hence the system 




Figure 7.3: The model results for the pionic HBT radii. The calculation uses the two- 
particle method and includes the Coulomb effects. The values of the model parameters 
are the same as in Fig. 17.11 

at freeze-out is, as expected, azimuthally symmetric in the transverse plane. We 
also note that the shape of the isotherms is consistent with the result presented in 
Fig. 4 of Ref. [80]. 

In Fig. 17.21 we give the model transverse-momentum spectra of hadrons. Com- 
pared to the RHIC results from Fig. 16.21 we find much larger multiplicities of the 
produced hadrons and smaller slopes of the spectra, indicating the larger transverse 
flow that is caused by the larger initial temperature. 

Our model calculations of the HBT radii are shown in Fig. I7.3[ The increase of the 
central temperature from Ti = 320 MeV to Ti = 400 — 500 MeV makes all the radii 
moderately larger. The ratio i?out/-Rside decreases by about 10% which is an effect 
of the larger transverse flow caused, in turn, by the larger initial temperature. 

7.2 Non-central collisions 

In this Section we present our results describing the peripheral collisions with Ti = 400, 
450 and 500 MeV, Tf = 150 MeV, and b = 7.6 fm. In Fig. Owe show the freeze-out 
curves. One can observe that the system is initially elongated along the y axis, but 
in the end of the evolution it becomes elongated along the x axis. This behavior is 
indicated by the crossing of the freeze-out curves. The change of shape is caused 
by the strong flow which transforms the initial "almond" into a "pumpkin" |llj . 
Azimuthally sensitive HBT probes such non-trivial behavior and can be used as a 
precise confirmation test for the existence of such effects in the data. 

In Fig. 17.51 we show the model transverse-momentum spectra of hadrons for the 
same values of the parameters. Compared to the RHIC case with Ti = 320 MeV and 
the same values of Tf, we find flatter spectra with higher multiplicity. In Fig. 17.61 
we show our results for the elliptic-flow coefficient. The stronger transverse flow 
generated in this case induces larger splitting between the pion-|-kaon V2 and the 
proton V2, with the values of the pion+kaon elliptic flow very similar to those found 
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Figure 7.4: The freeze-out curves for peripheral collisions at LHC, c = 20 — 30%, 
b = 7.6 fm, Ti = 1.0 fm, Ti = 400,450 and 500 MeV, and Tf = 150 MeV. The solid 
(dashed) line shows the in-plane (out-of-plane) profile. 
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Figure 7.5: The model results for the transverse-momentum spectra of tt~^ , K~^, and 
protons. The values of the model parameters are the same as in Fig. 17.41 

in the case T[ = 320 MeV. This result indicates the saturation of the elliptic flow of 
light particles for a given initial space asymmetry. On the other hand, the proton 
elliptic flow is significantly reduced. This observation is consistent with the findings 
of Kestin and Heinz discussed in Ref . [83] . Finally, in Fig. 17.71 we show our model 
calculations of the HBT radii. We note that the ratio i?out/-Rside is very close to one. 
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Figure 7.6: The elliptic flow coefficient V2- The parameters are the same as in Fig. 17.41 




Figure 7.7: The pionic HBT radii. The values of the model parameters are the same 
as in Fig. 17.41 



Chapter 8 



Uniform description of the RHIC 
data with Gaussian initial 
conditions 



In this Chapter we return to the discussion of the RHIC results obtained at the 
top beam energy of ^/sNN = 200 GeV. In the approach presented in Chapter E] we 
used the standard initial conditions and were able to reproduce the experimental 
data with the accuracy of about 10-20%. In the following we are going to show 
that further improvement in the description of the data may be achieved if one uses 
modified initial conditions, namely, we suggest to use a two-dimensional Gaussian 
as an initial profile for the energy density in the transverse plane. 

Of course, the elaboration of the precise initial conditions for hydrodynamics 
remains a challenge. Ideally, such initial conditions should be provided by the early- 
stage dynamics, for instance by the Color Glass Condensate [611 EHl El] ■ In practice, 
the theory of the early stage carries some uncertainty in its parameters. Moreover, 
different phenomena may be present in the early stage, see e.g. [HHIIHS], so its precise 
modeling is difficult. 



In the approach discussed in this Section we use a simple parameterization of the ini- 
tial energy density profile in the transverse plane. At the proper time of = 0.25 fm 
we assume that it may be described by a Gaussian 



where the values of the a and b width parameters depend on the centrality class 
and are obtained by matching to the results for (x^) and (y^) from GLISSANDO 
[57] , which implements the eccentricity fiuctuations of the system in the Glauber 
approach [HZl [HE]- The expression (a^ + 6^)/2 describes the overall transverse size 



8.1 Early start of hydrodynamics 




(8.1) 
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Table 8.1: Initial central temperatures and shape parameters for various centrality 
classes. 



c[%] 


- 5 


5 - 10 


10 - 20 


20 - 30 


30 - 40 


40 - 50 


T fMpVl 

J. I 1 iVlt; V J 


5nn 

o\jyj 


4Q0 


475 


460 


430 


390 


n TTTI 
Uj liii 




9 536 


2 284 


2 000 

^ .WWW 


1 771 


1 577 


U 1111 




2 84Q 


2 738 


2 5Q1 


2 448 

^ . 1 1(J 


2 305 


c[%] 


50 - 60 


60 - 70 


70 - 80 


80 - 90 


90 - 100 


30 - 80 


Ti [MeV] 


345 


303 


260 






330 


a [fm] 


1.401 


1.224 


1.040 


0.841 


0.511 


1.592 


b [fm] 


2.164 


2.019 


1.854 


1.649 


1.228 


2.307 



of the system, whereas (6^ — a^)/{h'^ + a?) parameterizes its eccentricity. The use 
of the Gaussian profile (18.11) leads to a faster development of the transverse flow as 
compared to the Glauber initial conditions. This effect is crucial for the success of 
reproducing the HBT radii. 

The values for centrality classes used in this work are collected in Table 18.11 
Since the equation of state is known, the initial central temperature Ti defines the 
initial central energy density ^i. It depends on centrality and has been adjusted 
to reproduce the total pion multiplicity. With the Gaussian initial conditions we 
continue the hydrodynamic evolution till the freeze-out temperature Tf = 145 MeV, 
where the system decouples and hadrons (stable and resonances) are generated ac- 
cording to the standard Cooper- Frye formalism [69] - see Chapters [6] and [71 Eq. 15.91 
The lower value of Tf leads to larger sizes of the system and larger transverse flow. 
This in turn results in flatter spectra and larger splitting of the V2 of pions/kaons 
and protons. 

The freeze-out hypersurfaces for two centrality classes are shown in Fig. 18.11 For 
comparison, in Figs. 16.1] and 16.41 in Chapter [6] we showed the hypersurfaces obtained 
in a calculation with the standard Glauber initial condition (and with same equation 
of state). We note that the Glauber initial condition results in hypersurfaces of a 
smaller transverse size and a longer evolution time, which translates into lower -Rgide 
and larger -Rout- More generally, we note a quite different shape of the hypersurfaces 
in Fig. 18.11 however, the flow values are similar, which results from the fact that the 
model parameters are adjusted in such a way that the slopes of the p^-spectra and 
the V2 coefficient are reproduced. This leads to practically the same values of the 
flow velocity at freeze-out for the two compared calculations. 
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Figure 8.1: Freeze-out hypersurfaces for c = — 5% (left) and c = 20 — 30% (right), ob- 
tained with the initial temperature T[ = 500 MeV and 460 MeV, for the central and 
semi-peripheral collisions, respectively, Ti = 0.25 fm and with the final temperature 
Tf = 145 MeV. The two curves show the in-plane and out-of-plane sections. The labels 
indicate the transverse flow velocity. 



We should also stress out the similarity of the volume-emission parts (with time- 
like normal vectors) to the blast- wave parameterization [89], with a short lifetime: 
about 9 fm for central (c = 0-5%) and 7 fm for mid-peripheral (c = 20-30%) colli- 
sions. However, the freeze-out hypersurfaces contain also the surface emission parts 
(with space-like normal vectors), absent in the traditional blast- wave parameteri- 
zations. We have checked that the potential problems from the non-causal surface 
emission [901 EH [92], are negligible, as less than 0.5% of particles are emitted back 
into the hydrodynamic region. The reason for this very small fraction is the sizable 
transverse flow velocity at large radii, as indicated by labels in Fig. 18. 11 which pushes 
the particles outward, as well as the fact that the hypersurfaces are not bent back 
at low values of t. We have determined that about half of the produced particles 
comes from the volume part and about half from the surface part of the freeze-out 
hypersurface. The surface emission is crucial for fitting the HBT data from RHIC, 
as also advocated in [221 E3] • 

In Fig. 18.21 we present our results for the hadron transverse-momentum spectra 
for positive pions and kaons and protons in the central and semi-peripheral colli- 
sions with the experimental data from PHENIX [77J. The bottom plot of Fig. 18.21 
represents the elliptic flow coefficient f 2 of joined pions -|- kaons and protons as a 
function of together with the PHENIX [7H] data points. We see a very good 
agreement of the one-particle observables. The V2 description of protons is better 
than the one presented in the Chapter El The splitting between the summed f 2 of 
pions and kaons is larger which is a sign of a larger flow difference that developed on 
the hypersurface. This description can be further improved by including the elastic 
rescattering processes among hadrons after freeze-out. 
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Figure 8.2: The transverse-momentum spectra of pions, kaons and protons for c=0-5% 
(left), c=20-30% (right), and the elHptic flow coefficient V2 for c=20-40% (bottom), plotted 
as functions of the transverse momentum and compared to the data from [771 178j . Model 
parameters same as in Fig. 18.11 

Figure 18.31 shows the HBT radii -Rside , -Rout , -Riong and the ratio -RoutZ-Rside- 
For the first time we can present a set of soft observables that are fully consistent 
with the experimental data. The most significant point is the correct description of 
the ratio i?out /-Rside in both central and non-central collisions. As in the previous 
calculation showed in Chapter [6], -Riong exceeds the data which is most likely caused 
by the boost-invariance of our model and can be fixed by implementing a rapidity 
dependent initial distribution. 

Although the success of the use of the Gaussian initial condition is evident, it is 
clear that a more detailed model describing the initial non-equilibrated stage of the 
evolution is necessary to address the issue of the microscopic validity of the initial 
Gaussian conditions. We leave this problem for our future investigations. 



8.2 Hydrodynamics preceded by free-streaming 

The very successful description of the soft observables that was presented in the pre- 
vious Section demanded a very early start of the hydrodynamic evolution, Xj = 0.25 
fm. Such short starting times have become now a common practice [OS]- Neverthe- 
less, the processes needed to thermalize the system require at least a few collisions 



8.2. HYDRODYNAMICS PRECEDED BY FREE-STREAMING 



65 




0.2 0.25 0.3 0.35 0.4 0.45 0.5 0.2 0.25 0.3 0.35 0.4 0.45 0.5 

kr [GeV] kj [GeV] 



Figure 8.3: The pion HBT radii i?sidc , Rout , -^long) and the ratio i?out/^sidc ^or central 
coUisions, compared to the data from [76]. Model parameters same as in Fig. 18.11 

to create local equilibrium. In this respect, the starting time of 0.25 fm seems to be 
very difficult to explain microscopically. 

In Chapter 14.31 we have introduced the process of parton free-streaming. The 
basic idea behind it is the approximation of the early, stage evolution of a partonic 
system formed in heavy-ion collisions. In our approach the parton free streaming 
is followed by sudden equilibration (described by the Landau matching conditions), 
and then by hydrodynamics. This kind of a formalism was proposed by Kolb, 
Sollfrank, and Heinz [3] several years ago in the context of the development of az- 
imuthally asymmetric flow. The approach assumes a sudden but delayed transition 
from a non-equilibrium initial state to a fully thermalized fluid. 

The process of free-streaming decreases the spatial asymmetry of the density 
profile. This effect was the reason for not implementing the free streaming before 
hydrodynamics. The less asymmetric initial condition for the hydrodynamic stage 
produces smaller amounts of the elliptic flow. Free streaming itself cannot generate 
azimuthal asymmetry in the momentum distribution. Interactions among produced 
particles are needed to achieve this goal. But a sudden equilibration preceded by 
free-streaming is in fact capable of developing azimuthally asymmetric flow. The 
energy-momentum tensor of the system changes abruptly into a diagonal form (in 
the reference frame co-moving with the fluid element). In this way the space-flow 
correlations are induced, which results in a collective elliptic flow, further enhanced 
by the subsequent hydrodynamic evolution. 

Our results presented in this Section were obtained with the free-streaming stage 
lasting from tq = 0.25 to Tj = 1.0 fm, followed by a sudden equilibration and 
hydrodynamics. The initial hydrodynamic time was shifted to a later time of 1.0 
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Figure 8.4: Freeze-out hypersurfaces for c = — 5% (left) and c = 20 — 30% (right). Cal- 
culations include the free-streaming phase from tq = 0.25 fm to Ti = 1.0 fm. The initial 
central temperature T[ equals 331 and 305 MeV, for central and non-central collisions, 
respectively. The freeze-out temperature is set to Tf = 145 MeV. 




Figure 8.5: The transverse-momentum spectra of pions, kaons and protons for c=0-5% 
(upper left panel), c=20-30% (upper right panel) and the elliptic flow coefficient V2 for 
c=20-40% (lower panel). Model parameters as in Fig. 18.41 Data from [771 178j . 
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fm. Such a change requires a different normahzation of the energy density at the 
center of the fireball, since during the time of 0.75 fm the system expands and its 
density decreases. To make an estimate of the initial energy at a later starting time 
we use the Bjorken scaling for the entropy |17]. We calculate the entropy density sq 
at To from the energy density Eq and scale it by the following formula 

■50^0 = SiTi = const. —>■ Si = sq — . (8.2) 

n 

As a result the initial central temperature for the most central collision (c = — 5%) 
is changed from T; = 500 MeV down to 331 MeV and for the peripheral case 
(c = 20 — 30%) the temperature is decreased from 460 MeV to 305 MeV. This type 
of scaling works surprisingly well as the method of construction the initial conditions 
at Ti = 1.0 fm. If we compare the hypersurfaces from the previous section, Fig. 18. H 
with the new ones obtained with the free-streaming. Fig. 18.41 then one can notice 
that they are almost identical. Both the shapes (starting from r = 1.0 fm) and 
values of the velocity on the hypersurface are very similar. 

The results of the calculations with a pre-hydrodynamics stage are shown in 
Figs. 18.51 and 18. 6[ As in Section 18.11 we have started from the Gaussian profile 
(18. ip . We notice very similar results to the ones obtained in the previous Section, 
not to mention the very good description of the data. Larger free-streaming times 
(r— To ~ 1.5 fm) spoil this agreement, as the flow becomes too strong. Again we have 
achieved a uniform agreement for soft physics at RHIC. In particular, the transverse- 
momentum spectra, the elliptic-flow, and the HBT correlation radii, including the 
notorious ratio -Rout/-Rside, are all properly described. 
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Figure 8.6: The pion HBT radii i?sidc , Rout , Riong, and the ratio i?out/-Rside for central 
colhsions, compared to the data from [76j. Model parameters same as in Fig. 18. 4[ 



Chapter 9 
Summary 



Our recently developed 2+1 (boost-invariant) hydrodynamic model has been pre- 
sented and used to i) describe the soft hadronic data collected in the central region 
of the relativistic heavy-ion collisions at RHIC and ii) to make predictions for the 
heavy-ion collisions at the LHC energies. We have addressed both the one- and two- 
particle observables: the transverse momentum spectra, the elliptic flow coefficient 
V2, and the pion HBT radii. The realistic equation of state for strongly interacting 
matter has been constructed that interpolates between the hadron gas model and 
the results of the QCD lattice simulations. The computational platform has been 
constructed, which combines the results of our hydrodynamic code with the sta- 
tistical hadronization model THERMINATOR. The satisfactory description of the soft 
hadronic RHIC data has been achieved with the standard initial conditions obtained 
from the optical limit of the Glauber model. Predictions for the future heavy-ion 
collisions at LHC have been formulated. The solution of the RHIC HBT puzzle has 
been proposed. It suggests the use of the modified Gaussian-type initial conditions 
for the energy density in the transverse plane. Finally, the processes of the free 
streaming of partons followed by the sudden equilibration have been incorporated 
in the model. The inclusion of the free-streaming stage allows for the delayed start 
of the hydrodynamic evolution, which is a desirable effect in the context of the early 
thermalization problem. 
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CHAPTER 9. SUMMARY 



Appendix A 

Relativistic thermodynamics of 
perfect gases 

A.l Grand canonical potential Q 

The first law of tliermodynamics for systems described by the grand canonical en- 
semble reads 

dE = TdS - PdV + fidN, (A.l) 

where E is the energy, 5* is the entropy, P stands for pressure, and N is the number 
of particles in the system (for relativistic systems the number of particles is not 
conserved and N should be interpreted as the baryon number or any other conserved 
quantity, for example strangeness). Energy is an extensive function of three extensive 
variables, E = E{S, V, N). By scaling the entropy, volume, and particle number by 
a factor A, we obtain the energy rescaled by the same factor 

E{XS, XV, XN) = XE{S, V, N). (A.2) 

By calculating the total derivative with respect to the scaling parameter A and using 
the relation ( lA.ip we find 



dE 
'dS 



v,N 9V 



N = E, 

sy 



(A.3) 



which leads to 

TS - PV + 12N = E. (A.4) 

The grand canonical potential Q is defined as 

n = E -TS - fiN = -PV. (A.5) 

To calculate the change of the potential Q we calculate the total derivative from Eq. 
(lA.Sp and use (lA.ip . In this way we find 



dn = -SdT - PdV - Ndfx. (A.6) 
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A. 1.1 The grand canonical partition function 

It is useful to start our consideration with the introduction of the partition function 
for perfect gases. In this case the grand canonical partition function is defined in 
the following way 

oo oo oo 

1=1 i=l ni=0 

where the index i runs over all microstates of the system and Zi is an individual 
partition function for that microstate. The index Ui is the number of particles 
with the energy in the state i. In the cases where the product is taken over the 
states with equal energy rather then over all microscopic states, we introduce the 
degeneration factor Qi giving the number of states with the same energy 



oo 



Z = l[ {zr ■ (A.8) 



i=l 



For indistinguishable particles satisfying the Bose-Einstein or Fermi-Dirac statistics 
we are able to evaluate explicitly the sum in Eq. (1A.7I1 



oo 



e T = __ (Bose-Einstein), (A. 9) 

ni=0 1 — e T 



1 / ^ 

e T = M + ) (Fermi-Dirac). (A. 10) 



ni=0 

A. 1.2 Quantum statistics 

From the statistical physics we know the connection between Z and Q, 

fi(T, V,/i) = -TlnZ, (A. 11) 

which gives 

oo oo 

n= -T\nZ = -Tin JJ(zi)^' = -Tj^di^'^^i- (A- 12) 

1=1 1=1 

This formula can be further rewritten in the form including Bose-Einstein (e = — 1) 
or Fermi-Dirac (e = +1) statistics, Eqs. (1A.9P and (lA.lOp . 

oo 

n = -Te ^^iln(l + ee^) . (A.13) 

i=l 

In the classical approach we can replace the sum in Eq. flA.lSp by the integration 
over the one-particle phase-space. This will lead us to the general form of the grand 
potential 

n{T,V,fi) = -TVge J -^ln(l + ee^). (A.14) 
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Here we assumed that all the degeneracies are equal Qi = g - this is the typical case 
where g is related to the spin degeneracy where g = 2s + 1. 

In the following we shall assume that the condition ix — E <{] is fulfilled. This 
condition is always true for bosons and it is also valid for fermions in the case n < m. 
In this case we use the series expansion of the logarithm, 

°° (—lY 

ln(l+x) = -^^ '-x^, (A.15) 

K=l 

to rewrite Eq. (1A.14I1 in the form 

n^siT,V,,)=TVge [ £ ^e^^ (A.16) 



A. 1.3 Boltzmann classical limit 

For particles which have small average occupation numbers one can use the approx- 
imation of the classical Boltzmann gas. This assumption is realized by using only 
the leading term in the expansion ( 1A.16I) . that corresponds to k = 1, 

n^UT,V,fi) = -TVg [ ^e^. (A.17) 

J (27r) 

We note that this approximation is independent of quantum statistics since in Eq. 
(1A.16I) we have always the term e'^'^^ which equals 1 for k = 1. 



A. 1.4 Massive particles 

The relativistic relation between the particle energy Ep, momentum p, and the rest 
mass m is £'p = \/p'^ + m^. Introducing this relation into Eq. ( 1A.16P enables us to 
calculate the potential Vt for massive particles 



Changing to the spherical coordinates and performing simple calculations we find 
^iT,V,^^) = ^^TVgeY,^^e^^ H dpp' e'^^^ (A.19) 



Let us now consider only the integral on the right-hand-side of Eq. (lA.lQp . By 
changing the integration variables from pio x with the substitution m?+p^ = m? x^, 
this integral takes the following form 

/oo I — 9 o /"OO -, 

dppe T = rrf j dxx[x -l)^e t ^. (A.20) 
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Table A.l: Grand canonical potential Q for various choices of m. 

quantum statistics 



m ^ 



m = 



■T'^VqeLL 



classical limit 









1 




(?) 






m = 









The formula above may be reduced to the integral definition of the modified Bessel 
function of the second kind Ky{z), see Eq. (ID .31) with z/ = 2 and z = m k/T, 



m m K 



3 T 



ctx — Ij ^ e T 



2 



(A.21) 



Thus we write the final form of the grand canonical potential for massive particles 
with finite chemical potential and quantum statistics as follows 



Q{T,V,fi) 



271- 



■T^Vm'^geY^ 

K=l 



'^e^" K2 [j^ ) (m > 0, /i < m). (A.22) 



A. 1.5 Massless particles 

Of the special interest is the case of massless particles. Such a limiting case is 
appropriate for gluons and also for quarks at sufficiently high temperatures 0. The 
thermodynamical properties of massless particles follow directly from the formula 
(1A.14P where we set m = and /x = 0. The alternative way is to consider the 
limiting case of Eq. (1A.22I) . For massless particles the modified Bessel function 
Ki,{z) in Eq. ( ]A.22I) may be expanded into a power series in the vicinity of the 
point m = according to Eqs. flD.7l) . This leads us to the expression 



niT,V,,) = ^^T^Vm^geJ2^-^e^ 



K = l 



(A.23) 



^Typically the light quarks u and d may be regarded as massless particles. If the temperature 
is sufficiently high, in addition the strange quark may be regarded as a massless particles. 
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Table A. 2: Mean number of particles N for various choices of m. 

quantum statistics 



TO 7^ 



TO = 



TO 7^ 
TO = 



K=l 



classical limit 



— -TV qm^ eT K2[ — 



Since only the first term of this expansion gives non zero contribution in the limit 
m — > 0, the grand potential for massless particles reads 



n{T,V,^^) = -T'VgeY, 



^T^VgeLu (-eeT 



(A.24) 



Here we have introduced the polylogarithm function defined by the series 

CO ^ 



k=l 

In the limit of vanishing chemical potential the values Li4(l) (required for bosons) 
and Li4(+1) (required for fermions) should be obtained from the analytic continua- 
tion which gives Li4(l) = 7r^/90 and Li4(— 1) = — 77r^/720. In this way we obtain 

n{T,V,fi = 0,m = 0) 



and 



Q{T,V,fi = 0,m = 0) 



90 
77r2 



720 



gVT^ 
gVT' 



(Bose-Einstein) (A.26) 
(Fermi-Dirac). (A. 27) 



A. 1.6 Vanishing chemical potential 

Sometimes one assumes that the chemical potential is very small and may be ne- 
glected. In this case the grand canonical potential, Eq. (lA.16p . takes the form 



n{T,V,fi = 0) = -TVge 



E 



-er _Ep 



(A.28) 
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or in the case of classical massless gas with no chemical potential 

nnL{T,V,fi = 0) = -^gVT\ 

A. 2 Other thermodynamic variables 



(A.29) 



Given the formula for the differential of the grand potential flA.6p we can evaluate 
the mean number of particles using the following expression 



N 



9/i 



(A.30) 



With the most general definition of Q, valid for massive gas with quantum statistics, 
see Eq. (1X221) . we find 



NiT,V,fx) 



27r' 



T V vr? g e 



K = l 



<-^'%*-/<-Jf.). (A,3i; 



Neglecting the quantum statistics will lead to the Boltzmann classical approximation 
in the form 

ATblIT, V,fi) = ^TVe^m'gK2 (^) . (A.32) 

In the further step we can do the massless-gas assumption and following the same 
formalism as in subsection IA.1.51 one gets 



Nj,^{T,V,ix) = -T'gVe^. 



The formula for the entropy follows from the equation 



which for massive particles obeying the quantum statistics gives 
1 



—e 



K=l 



Similarly, to calculate the pressure we use the formula 

P = 



dV 



and get 



K = l 



(A.33) 



(A.34) 



{AT-^)K2(-K^+mKK^(^-K] . (A.35) 



(A.36) 



(A.37) 



A.2. OTHER THERMODYNAMIC VARIABLES 



Table A. 3: Entropy for various choices of m. 

quantum statistics 



m ^ 



m = 



classical limit 



m ^ 
m = 



Table A. 4: Pressure for various choices of m. 

quantum statistics 







m = 



1 

'2^ 



T'^gm'^e 



K=l 



T 



classical limit 



m ^ 
m = 



^^T^gm^e^K,{ 



t) 



Finally, the energy is obtained from the relation 

E^TS-PV + iiN 

which gives 



(A 



3TK2 [-Kj +mKKi i^-n 



(A 
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Table A. 5: Energy for various choices of m. 

quantum statistics 



TO ^ 



TO = 



«=1 



-^T^VgeLu[-ee^) 
classical limit 



TO ^ 

TO = 



27r2 



TVgm^eT 



Appendix B 

Properties of relativistic 
hydrodynamic equations for /i^ = 



This Appendix presents the details of transformations of the relativistic hydrody- 
namic equations which were applied to obtain the final form used in the numerical 
calculations. It also introduces the basic hydrodynamic definitions that are used in 
our work. 



B.l Basic definitions 

The fluid four-velocity is defined as 



= (ii^iiO =7(l,v) ^-f{l,v^,Vy,V:,), (B.l) 
where 7 is the Lorentz factor defined as 

The derivative of 7 equals 

^7 = -f^v'dv' = -fhdv, (B.3) 

where 

vdv = ^~^d'y, V = '^v^ + v^ + v'^ = ^v^^-\-v^. (B.4) 

Here v± is the magnitude of the transverse velocity vj^. The total time derivative 
in the Cartesian coordinates reads 

d d , d d d d d /t-, ^\ 

while in the cylindrical coordinates it has the form 

d d ^ d ^v±sma d ^ d (B 6) 
dt dt ^ dr r d(f) dz 

Here a is the angle between the transverse velocity and radial velocity v^, 
Vr = v± COS a, v^ = v± sin a. 
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B.2 Covariant form 

The conservation of the energy-momentum tensor of the perfect fluid in its most 
general form yields 

[{e + P) u^u" - P g''"] = 0. (B.7) 

We transform Eq. ( IB .71) using thermodynamic identities e + P = T s and dP = s dT 

to get 

Tu'' + SM^ d^iTu"") = sd^T. (B.8) 

Multiplying both sides of Eq. (IB. 80 by Ui, and using the normalization condition for 
the four- velocity, u^Uiy = 1, we arrive at the following formula 

Tdf,{su^) + su^Tu^d^ + su^d^^T = su^d^T. (B.9) 

The term that includes Uud^u" is equal to zero because of the four-velocity normal- 
ization condition. Thus, we find the expression 

Td^isu^") + su^'d^T = su^d^T. (B.IO) 

This equation directly implies that 

d^isu'') = 0. (B.ll) 

Eq. (IB. lip states that the entropy is conserved in the system - the hydrodynamic 
expansion is adiabatic. 

In the next step we return to Eq. ( 1B.8P and use Eq. ( IB.lip to rewrite it in the 
following form 

u^'d^iTu") = d^T. (B.12) 

Eq. ( IB. 121) is the acceleration equation - a relativistic generalization of the Euler 
equation known from the classical fluid dynamics. Equation ( IB. 121) includes only tree 
non-trivial equations. It is easy to check this property if we project this equation on 
the four- velocity u^. The final form of the hydrodynamic equations in the covariant 
form consists of a set of four equations 

u'^df.iTu'') = d^T, 

d^isu^) = 0, (B.13) 

which should be supplemented by the equation of state. 



B.3 Non- covariant form 

In this Section we show how to rewrite relativistic hydrodynamic equations in the 
non-covariant form. We start with the acceleration equation in the covariant form, 
Eq. ( IB. 121) . and show that the corresponding non-covariant form is 



— (T7v) + V(T7) = V X [V X (T7V)] 



(B.14) 
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In order to prove this property we write the ith (spatial) component of Eq. flB.141) 
as 



d_ 

di 



(B.15) 



The next step is to contract the e tensors with the help of the formula 

and to carry out the sums with the Kronecker delta. In this way we find 

d 

Q^iT-fv') + V*(T7) = v^V'{T-fv^) - v^V^{T-fv'). 



(B.16) 



We move the terms in Eq. ( ]B.16p that include T'yv^ to the left-hand-side and 
multiply both sides of this equation by the Lorentz 7 factor. This procedure yields 
the equation whose sides are given by the expressions 



LHS = 7 ^{T-fv') +-fv^V^{T-fv') = u'^doiTu') + u^djiTu') = u''d^{Tu'), 



RHS = -iv^V\T-iv^)--iV\T-i)=-i^{v^ -l)V'T + T 
= -V'T = d'T. 



7(t;2 



+ i 

7. 



V'7 



One observes that both sides of Eq. (IB. 120 are reproduced after those manipulations. 
The non-covariant form of the entropy conservation law, Eq. (IB.lip . may be written 
in the straightforward way as 



d_ 
di 



(s7) + V(s7v) = 0. 



(B.17) 



B.4 Temperature equation 



In this section we study in more detail the temperature equation flB.141) . We start 
with the iih component of Eq. flB.14p . which in the tensor notation has the form 



d_ 
di 



{T-fv') + V\T-f) = v^V'{T-fv^) - v^V\T-fv'). 



(B.18) 



We expand separately the left- and the right-hand-side of Eq. (1B.18P to find 
LHS 



,dT ^ .97 _ dv' dT ^d-f 



dt 
.dT 



dt dt dx^ dx 



.dv^ 



+ + 7^ + T-i^v^ — 

dt dx' dx^ 



T7 



■ d\nT 
~dt 



+ 7 t; f 



dt ^ dt 



+ 



2 „,fc 



dx^ 



+ 71; 



dx^ 



[B.19) 
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RHS 



■ dT ■ ■ d'y dv^ ■ dT ■ ■ d'j dv^ 

v^'yv^ —— + v^Tv^ —— + v^T'y—— — v^'yv'^-;—- — v^Tv^—— — v^T^y—— 
ax* ax* ax* ox^ ox^ ox^ 

(9x* (9x* 9x* dx^ ' dx^ ' dx^ 

V 



^ ^dT ^ 3 J i kdv'' „ jdv' 
I 7 t>-'t> V — — r — I 'yv' 



T7 



,d\nT 



9x* 9x* dx 

By comparing LHS with RHS we get 
d\nT 



dx^ 



— 



dx^ 



(9x* 



d\nT ,(91nT 
h 



dt 



dx^ 



dv"- ,(9f* 2 k 



ht;^ 



dt 



dx^ 



The final form of the temperature equations in the Cartesian coordinates is 

SlnT 



^ ' dx' 



■ dlnT dv^ 2 i kdv^ q 
dt dt dt ^ 



(B.20) 



or 



^SlnT ^ dlnT ^ dvx ^ 2 dv 

dx dt dt dt 

9,(91nT d\n.T dvy <, dv 

2 (91nT dlnT dv^ 2 dv 



0, 
0, 
0. 



(B.21) 
[B.22) 

[B.23) 



For further use it is convenient to consider linear combinations of Eqs. flB.2ip - 
flB.23p . At first we multiply Eqs. flB.2ip - flB.23p by the appropriate components 
of the velocity (for example, Eq. ( IB. 211) is multiplied by Vx, etc.) and add them 
together. The second combination is obtained if we multiply Eq. ( IB. 211) by Vy and 
subtract Eq. (]B.22|) multiplied by Vx- As the third independent equation we take 
the unchanged formula (]B.23|) . In this way one finds 



d\nT 



11 



91nT rflnT 

+ 



dx 



dt dt 
dlnT 

Vx^ 

dy 



dlnT dlnT dv. 



dz 



dt ~'~ dt 



1 dv 




(B.24) 


^ l-v^'^dt 


0, 


dvx dvy 




(B.25) 


0, 


Vz dv 




(B.26) 


^ 1-v^^ dt 


0. 
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These equations can be also rewritten in the cyhndrical coordinates in the form 

, n . dlnT dlnT 1 dv 



(B.27) 



dlnT VrCosa d\iCLT\ ^ f da frSina' 

,l-v)\VrSma— H3— -ttH 1 = 0' 

or r om I \ dt r 



(B.28) 

9, (91nT c/lnT dv^ Vz dv 

(B.29) 

We note that the total derivative in the above equations has the form of Eq. (lB.6p . 



B.4.1 Boost-invariance 

For the boost invariant systems the longitudinal velocity has the form Vz = z/t. We 
may insert this form into the hydrodynamic equations, calculate all the necessary 
derivatives and then set z = 0. In this way we obtain the hydrodynamic equations 
in the plane z = 0. 

,9 .dlnT dlnT 1 dvr ,^ 
n, / 91nT 91nT\ dv^ dvy 

Again we can make use of the cylindrical coordinates to get 



ainT dlnT 1 dv 



'1 



9 / dlnT t;,. cosa (91nT\ ^ ( da t;^. sina 

\-V,)\Vr'mia — — — - -TT H I = 0- 

or r Off) \ dt r 



(B.32) 



fB.33) 



The total derivatives for the boost-invariant systems are reduced at ^ = to the 
expressions 

It = I + "-^ + "4 (Cartesian), (B.34) 
d d d V sin a d 

— = — + VrCosa— — I — — (cylindrical). (B.35) 

dt dt or r o<p 
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B.4.2 Cylindrical symmetry 



In certain physical situations the system described by the relativistic hydrodynamics 
has cyhndrical symmetry. In this case the hydrodynamic equations are reduced to 
the form 



iv'' - 1 



d\nT d\nT 



+ 



d\nT 

dz 



+ V. 



dt dt 
dlnT dv.. 



1 dv 



dt 



+ 



dt 



+ 



dv 



0. 



(B.36) 
(B.37) 



These equations have the same form in cylindrical coordinates. Combining both the 
boost-invariance and cylindrical symmetry gives us only one independent equation 



d\nT d\nT 1 dvr Vr dvr 

Vr—::— + — ^ + " ^"T^ + , . ^ = 0. 



dt 



dr 



1-v^ dt 



1 — v"^ dr 



(B.38) 



B.5 Entropy equation 



In this Section we discuss the entropy conservation law, Eq. ( IB. 171) . Following the 
same steps as in the previous Section we obtain 



ds d'-y ■ ds ■ d'j dv 



7S 



dins 



+ V 



dt dt 
:d\ns Ok 



dx^ 



+ v' 



dx 
dv'^ 



+ 



dx^ 
dv^ 
dx^ 



dt dx' V 9t dx' 

The final form of the entropy equation has the form 

dins V dv 9f* 
dt ^ 1 - v^~(U ^ dx' ^ ' 
which in the Cartesian and cylindrical coordinates has the form 



dins 
~dt 



+ 



dv dvr dv„ 



+ 



+ 



1 — v'^ dt dx dy 



+ 



dv, 
dz 



0, 



and 



c? In s V dv dvr sin a dvr 

dt 1 — f 2 dt dr r dS 



da Vr cos a f da 

Vr sm a— H ^ + ^ 

dr r \d(p 



dvz 
dz 



0, 
0. 



(B.39) 



(B.40) 



(B.41) 



B.5.1 Boost-invariance 

The z component of velocity equals = z/t, hence its derivative with the respect 
to that coordinate gives a term j. With this remark in mind we perform the rest of 
the calculations at 2; = and find 



dins Vr dvr dvr dv,, 1 

\ H H + - = 0, 

dt 1 — v'^ dt dx dy t 



(B.42) 
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or 

dins 1 dvr dvr 2 . da vl cos a f da \ Vr 

"^^ + 1^^- 9^-"^"""^;^ + ^— J^T = °- 



(B.43) 



B.5.2 Cylindrical symmetry 



For cylindrically symmetric systems entropy density and transverse velocities are 
independent of direction in transverse plane. The entropy equation with cylindrical 
symmetry has the form 

dins Vr dVj. 1 dVy Vr dvz „ , ,v 

at 1 — at 1 — or r oz 
where we have used the definition for the transverse velocity divergence 

dVx dVy dVr Vr 

^ + ^ = + -■ B.45 
ox Oy Or r 

When both the boost-invariance and cylindrical symmetry is applied to the entropy 
equation we get 

9 In S 9 In S Vr dVr 1 dVr Vr 1 ^ ,g 

dt ^ dr 1 — v^ dt 1 — v^ dr r t 
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Appendix C 



Notation 



In this Appendix we collect the symbols used in the Thesis. 



Thermodynamics 



grand canonical potential, 

Z or z grand canonical partition function, 

T temperature, 

V volume, 

E energy, 

E , . 

£ = — energy density, 

S entropy, 

S , . 

s = — entropy density, 

P pressure, 

N particle number, 

/X chemical potential, 

/iB baryon chemical potential, 

/is strange chemical potential, 

/ii3 isospin chemical potential, 

e statistics identifier (+1 for FD, -1 for BE), 

g degeneration factor, 

Cg sound velocity, 

Tc critical temperature, 

m„ pion mass. 
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APPENDIX a NOTATION 



Hydrodynamics 



T'^'^ energy momentum tensor, 

u** = 7 (1, v) four- velocity, 

= {vx,Vy,Vz) velocity, 

radial velocity, 
transverse velocity, 

a = ta.n~'^{vy/v.,.) — (f> angle between transverse velocity v_l and the radial velocity v^, 
r]± — tanh^^ v± transverse fluid rapidity, 
7 = (1 — Lorentz factor, 

baryon four-current, 
T proper time, 

$ thermodynamic potential used in the Baym formalism, 

j4+, A-,A auxiliary functions used in the Baym formalism. 



Initial conditions for hydrodynamics 



nuclear profile, 

parameters for Woods-Saxon nuclear profile, 

wounded nucleon profile, 

binary collisions profile, 

thickness function, 

total inelastic pp cross-section, 

total inelastic cross-section in Au Au, 

mixing factor between WN and BC nuclear profiles, 

impact vector, 

initial proper time, 

initial entropy density at r = 0, 

initial energy density at r = 0, 

initial temperature at r = 0, 

initial Hubble fiow. 



Freeze-out 



spatial azimuthal angle, 
angle in the t — p plane, 

spacetime rapidity, 

distance from the point {t = x = 0, y = 0) to the hypersurface 

point with coordinates (0, Q , 

transverse velocity on the hypersurface, 

angle between v± and Vr on the hypersurface, 

particle's transverse momentum, 

particle's azimuthal transverse momentum angle, 

particle's rapidity, 

particle's mass, 

three-dimensional element of the hypersurface. 



P 

Po,ro,a 

PWN 

Pbc 
Ta 

_AuAu 
K 

b 

Si 
Ei 

Ti 
Ho 



= W- 

2 t 
rf(<^,C) 



V_L 

a 

Pt or pj_ 

y 

m 



Appendix D 

Mathematical supplement 



D.l Modified Bessel functions 

The second-order ordinary differential equation 



^ , , +z^-{z' + u')w{z)^0 (D.l) 
az'^ az 

has a solution in form of the modified Bessel function of the first kind and 
second kind K,y{z) 

w(z) = cj,(z) + C2K,{z). (D.2) 

Modified Bessel function K^{z) has a simple integral representations through the 
exponential function and power functions in the integrand 

i^n(^)- ^.rV + |) yi Re(^) > -- A Re(z) > 0. (D.3) 

The Bessel function Ki,{z) satisfies the following recurrence identities 



(D.4) 



The derivative of the Bessel function Ky{z) has a rather simple and symmetrical 
representation that can be expressed through other Bessel Ky[z) functions with 
different indices 

= -\ (Ku-iiz) + K^+iiz)) , (D.5) 
this derivative can also be represented in other form 

= -K^-iiz) - -K,{z) = -K{z) - K+,{z). (D.6) 
oz z z 
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APPENDIX D. MATHEMATICAL SUPPLEMENT 



Series representations of the modified Bessel function Ki,{z) about the point z = 
is very complex. We will limit ourselves to those series expansions useful in our 
calculations 

K,{z) = - + 0{z'), 



(D.7) 



Mz)^^-l + 0{z% 

The asymptotic expansion of the modified Bessel function K^{z) about 2; — > oo has 
the form 

K(z)^e-^(J^z-'^+0[z-'^)] (D.8) 
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